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CHAPTER  1 


INTRODUCTION 

The  problem  of  predicting  electromagnetic  scattering  from  periodic  structures  has 
been  studied  for  well  over  a  century.  As  early  as  1786,  David  Rittenhouse,  an  American 
physicist,  published  his  experiments  on  the  diffraction  of  light  through  a  periodic  grating 
[1,  2],  Fraunhofer  constructed  some  of  the  first  practical  diffraction  gratings  in  1823  to 
split  a  beam  of  light  into  its  various  colors.  At  lower  frequencies.  Hertz  utihzed  a  metallic 
grid  in  his  experiments  in  1889.  Within  a  few  years,  several  theoretical  papers  had 
appeared  attempting  to  explain  Hertz’s  experiments,  the  first  by  J.  J.  Thompson  in  1893 
[3].  From  that  time  to  the  present,  the  problem  has  received  intermittent  attention,  driven 
primarily  by  the  discovery  of  new  applications,  the  development  of  printed  circuit 
techniques,  and  the  advent  of  the  high-speed  digital  computer.  The  last  of  these  factors  has 
permitted  a  wide  range  of  complicated  geometries  to  be  treated  efficiently  and  accurately 
using  a  handful  of  systematic  techniques. 

At  microwave  frequencies,  the  two-dimensional  versions  of  these  surfaces  are 
generated  by  replicating  a  conducting  patch  or  aperture  in  a  periodic  fashion  along  two 
surface  coordinates,  which  may  or  may  not  be  orthogonal.  The  shape  of  the  generating 
patch  or  aperture  is  used,  to  some  extent,  to  control  the  frequency  response  of  the  surface, 
for  example,  to  increase  the  bandwidth  or  to  make  the  frequency  response  less  sensitive  to 
incident  angle.  Special  considerations  apply  to  circular  polarization;  otherwise,  the  element 
shape  is  arbitrary. 

As  a  function  of  frequency,  electromagnetic  scattering  from  periodic  surfaces  is 
characterized  by  an  infinite  number  of  alternating  reflection  and  transmission  bands.  At 
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cenain  frequencies,  a  penodic  array  of  printed  elements  can  reflect  an  incident  field  as 
though  it  were  a  solid  perfectly-conducting  sheet.  At  other  frequencies,  the  surface  is 
perfectly  transparent.  It  is  this  property  that  is  most  often  used  in  satellite  antenna 
applications,  where  a  periodic  surface  is  used  to  generate  a  second  focus  for  a  parabolic 
reflector.  It  is  also  the  property  that  gives  it  the  name  frequency  selective  surface  (FSS). 

Nearly  all  of  the  theoretical  analyses  of  frequency  selective  surfaces  has  been  based 
on  the  assumption  that  the  surface  is  infinite  in  extent  along  both  coordinate  axes.  For 
plane-wave  incidence,  the  cell  currents  differ  only  by  the  known  progressive  phase  shift  of 
the  incident  field;  thus,  the  fundamental  unknown  is  either  the  current  density  or  aperture 
field  on  a  single  unit  cell.  When  frequency  selective  surfaces  are  implemented  in  real 
systems,  they  are  finite,  and  due  to  their  different  environments,  currents  flowing  in  the 
interior  regions  of  the  surface  can  be  quite  different  from  currents  flowing  near  the  edges. 
The  goal  of  this  research  is  to  study  the  effects  of  truncating  a  periodic  surface  to  a  finite 
size. 

A  frequency  selective  surface  that  has  been  truncated  to  finite  dimensions  is  not 
periodic.  Unless  otherwise  indicated  by  the  context,  the  term  “periodic  surface”  will 
always  be  used  to  refer  to  the  infinite  planar  surface,  which  is  (infinitely^  periodic  in  one 
dimension  (in  the  case  of  strip  gratings)  or  two  dimensions  (in  the  case  of  true  two- 
dimensional  elements).  The  term  “dichroic  surface”  will  be  avoided,  since  many 
applications  are  concerned  with  more  than  two  frequencies. 

This  thesis  is  written  from  the  problem- solver's  point  of  view.  Most  of  the  emphasis 
IS  on  the  formulation  of  the  equations  and  on  the  con.siderations  that  are  essential  to 
efficiently  and  accurately  solve  those  equations  on  the  computer.  Less  emphasis  will  be 
placed  on  physical  interpretation  of  the  results.  This  thesis  is  not  intended  to  be  a  survey  of 
the  theory  of  periodic  arrays.  A  complete  discussion  of  the  properties  of  the  Floquet  space 
harmonics,  surfaces  that  are  periodic  in  other  separable  coordinate  systems,  the  application 
of  Babinet's  principle  to  complementary  screens,  iterative  solutions  of  the  scattering 
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equations  via  the  conjugate  gradient  method,  and  scattering  matrices  for  cascading  arbitrary 
dielectric  layers  and  layers  of  printed  elements  are  all  treated  elsewhere  and  will  not  be 
repeated  here  [2,  4],  The  scope  of  this  thesis  includes  only  material  that  directly  pertains  to 
the  analysis  of  finite  frequency  selective  surfaces,  or  material  that  is  not  readily  available 
elsewhere. 

The  thesis  is  organized  as  follows:  Chapter  2  introduces  the  geometries  considered  in 
the  body  of  the  text,  identifies  the  fundamental  unknown  in  each  case,  and  summarizes  the 
relationships  between  the  fundamental  unknowns  and  the  total  current  density  in  the  plane 
of  the  printed  surface.  Chapter  3  provides  a  review  of  the  spectral-Galerkin  technique, 
with  a  detailed  look  at  the  derivation  of  the  spectral-domain  Green's  function  for  an 
arbitrary  /V-layer  dielectric  medium.  Chapter  4  highlights  some  of  the  theoretical  and 
numerical  considerations  that  are  essential  to  implementing  the  spectral-Galerkin  equation 
on  the  computer.  A  few  representative  numerical  results  are  presented  in  Chapter  5,  the 
response  of  frequency  selective  surfaces  to  arbitrary  incident  fields,  such  as  the  near  field 
of  an  electromagnetic  horn,  is  the  subject  of  Chapter  6,  and  Chapter  7  summarizes  the  main 
conclusions  and  suggests  directions  for  future  work. 
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CHAPTER  2 

REPRESENTATION  OF  FUNCTIONS  IN  PERIODIC  AND  APERIODIC 

DOMAINS 

A  few  years  ago,  the  spectral-domain  formulation  of  the  method  of  moments  was 
developed  and  successfully  applied  to  a  number  of  different  problems.  Although  the 
connection  is  not  always  clear,  many  of  the  formulations  for  predicting  electromagnetic 
scattering  from  periodic  surfaces  that  directly  make  use  of  the  fact  that  the  scattered  field  is 
composed  of  a  discrete  spectrum  of  plane  waves  are  essentially  equivalent  to  the  method 
presented  here.  Some  of  the  foundational  papers,  for  example,  Montgomery  [5],  Chen  [6], 
Chen  [7],  and  Tsao  and  Mittra  [8],  have  very  different  points  of  view,  terminology,  and 
notation,  but  they  all  arrive  at  the  same  fundamental  equation,  and  solve  it  using  the  same 
techniques.  In  the  next  chapter,  we  will  take  a  look  at  the  spectral-Galerkin  technique  as  it 
applies  to  a  fairly  broad  class  of  FSS-type  problems,  but  first,  in  order  to  avoid  repetition, 
it  is  necessary  to  review  some  basic  concepts  on  representing  functions  in  periodic  and 
aperiodic  domains 

One  of  the  most  useful  means  of  evaluating  the  effects  of  truncating  a  periodic  surface 
is  to  compare  the  currents  and  scattered  fields  obtained  for  the  finite  problem  with  the 
corresponding  results  for  the  periodic  problem.  The  amount  that  the  two  results  differ  is 
entirely  due  to  the  finite  versus  infinite  character  of  the  two  problems.  However,  in  order 
to  assign  all  of  the  difference  to  the  truncation  effect,  it  is  necessary  to  solve  both  problems, 
the  finite  problem  and  the  periodic  problem,  using  similar  approximations,  particularly  if 
one  is  interested  in  comparing,  for  example,  the  induced  currents  on  the  two  structures. 
For  this  reason,  both  finite  and  periodic  surfaces  will  be  considered,  concentrating  more 
heavily  on  the  finite-FSS  problem. 
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Having  stated  that  the  subject  of  the  thesis  is  electromagnetic  scattering  from  finite 
frequency  selective  surfaces,  the  subject  matter  is  still  too  broad  and  varied  to  allow  a 
comprehensive  treatment.  The  problems  considered  here  have  been  restricted  by  the 
following  simplifying  assumptions.  First,  panially  due  to  the  fact  that  the  method  of 
solution  is  a  Fourier-transform  technique,  only  planar  geometries  are  considered.  The 
problem  of  predicting  the  scattered  fields  from  a  frequency  selective  surface  which  has  true 
two-dimensional  curvature,  such  as  a  parabolic  dish  or  hyperbolic  subreflector,  is  still  a 
challenging  problem  computationally,  and  one  which  has  only  been  modeled  using  crude 
approximations  or  by  relying  heavily  on  experimental  data.  Cwik  has  discussed  surfaces 
which  are  periodic  in  cylindrical  and  spherical  coordinate  systems  [2],  and  he  has  also 
considered  finite  and  periodic  strip  gratings  on  the  surfaces  of  infinite  cylinders  [9]. 
However,  the  general  problem,  for  the  most  part,  has  not  been  addressed  -  nor  will  it  be 
addressed  here. 

Second,  although  the  method  is  not  limited  to  single  layers  of  metallization,  for  the 
purposes  of  this  thesis,  the  geometry  considered  is  a  single  layer  of  printed  elements 
embedded  in  an  A^-layer  dielectric  medium,  where  each  layer  is  assumed  to  be  uniform, 
homogeneous,  and  characterized  by  the  complex-valued  relative  dielectric  permittivity  £>. 
Under  these  assumptions,  the  problems  of  interest  here  fall  into  four  or  five  general 
categories,  which  are  treated  in  order  of  their  computational  difficulty. 

As  a  corollary  to  the  uniqueness  theorem,  the  fundamental  unknown  in  the  integral 
equation  can  be  chosen  as  either  the  induced  electnc  current  density  on  the  metallized 
portions  of  the  surface  or  as  the  tangential  electric  field  in  the  apertures.  In  the  periodic 
problem,  the  choice  is  somewhat  arbitrary,  but  it  is  usually  computationally  more 
convenient  to  choose  as  the  unknown  that  quantity  that  occupies  the  least  physical  area.  In 
the  finite- a  .y  problem,  the  choice  is  more  clear.  In  the  case  in  which  a  finite  array  of 
apenure?  . ,  c.>sumed  to  exist  in  an  infinite  conducting  sheet,  the  idea  of  representing  the 
current  deuii..  over  the  infinite  conducting  sheet  does  not  seem  to  be  an  option.  In  the 


6 


same  way,  for  a  surface  formed  by  printing  a  finite  number  of  conducting  patches  on  an 
infinite  dielectric  substrate,  choosing  the  tangential  electric  field  as  the  fundamental 
unknown  would  seem  to  be  impractical,  since  this  tot)  would  involve  representing  the 
unknown  over  an  infinite  area.  All  of  the  finite-FSS  examples  that  we  will  consider  are  of 
the  second  kind,  where  the  metallized  region  is  assumed  to  be  finite  in  some  sense,  and 
hence,  the  induced  current  density  on  the  surface  will  always  be  taken  as  the  fundamental 
unknown. 


2ol  Periodic  Strip  Grating 

The  first  two  types  of  problems  that  are  considered  are  strip  gratings  [9- 12];  infinitely 
long  layers  of  metallization  printed  on  a  dielectric  substrate,  as  shown  in  Fig.  2.1. 
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These  are  one-dimensional  structures  in  the  sense  that  the  geometry  is  independent  of  one 
of  the  surface  coordinates.  As  a  result,  they  are  among  the  simplest  structures  of  any 
practical  interest. 

In  the  fu^t  case,  we  will  consider  a  periodic  grating,  in  which  a  unit  strip  is  replicated 
in  a  periodic  fashion  an  infmite  number  of  times  in  the  direction  perpendicular  to  the  strips. 
The  metallized  areas  are  assumed  to  be  zero-thickness,  perfectly-conducting  regions  of 
width  a  and  periodicity  b,  embedded  in  N  uniform  homogeneous  dielectric  layers, 
described  by  complex  relative  permittivities  f>,  and  thicknesses  t,-. 

In  order  to  determine  how  this  structure  scatters  an  incident  electromagnetic  field,  the 
objective  is  to  first  write  an  expression  for  the  induced  current  density  on  the  entire  surface 
in  terms  of  the  current  on  a  single  element.  If  the  incident  field  is  assumed  to  be  a  plane 
wave  incident  from  the  direction  derived  from  the  magnetic  and  electric  vector 

potentials  A  and  F,  where 

F  =  (2.1) 

and, 

=  ^0  sin0j  COS0, 

ky  =  ^0  sin0,  sim^,  (2.2) 

k^  =  kQ  cosd; , 

then  the  y  variation  of  the  induced  surface  current  will  have  the  same  exponential  phase 
dependence  as  the  incident  field.  Adopting  the  notation  of  Papoulis  [13],  the  .x  and  y 
components  of  the  total  induced  current  density,  Jx  (Jt,y)  and  Jy  ix,y),  can  be  written  in  the 
form. 


J^ix,y)  = 
Jy(x,y)  = 


^  5{x  -  mb) 

m  =  -oo 

oo  j 

^  Six  -  mb)  r-  ^  V  , 


m  =  -<» 


(2.3) 
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or  in  vector  notation. 


JU,}')  =  j 


eJi^yy 


(2.4) 


With  the  y  dependence  removed,  the  remainder  of  the  total  current  density,  in  draces,  is  a 
function  of  x  only.  It  can  therefore  be  represented  as  a  one-dimensional  convolution  of 
Jo(x),  the  current  density  on  the  center  strip  along  the  cut  y  =  0,  and  an  infinite,  periodic 
sequence  of  delta  functions,  which  are  weighted  by  the  progressive  phase  shift  of  the 
incident  field. 

The  representation  given  by  Eqs.  (2.3)  and  (2.4)  is  not  unique.  Another  possibility 
would  have  been  to  factor  out  the  phase  shift  due  to  the  incident  field  in  both  x  and  y, 
leaving  a  periodic  function  of  x, 


J(jc,y)  =  \ 


JoOc)  *  ^ 


m  =  -<» 


(2.5) 


The  difference  between  the  two  representations  is  usually  insignificant,  but  it  is  possible  to 
imagine  situations  in  which  one  representation  might  offer  some  advantages  over  the  other. 
Anticipating  the  need  for  the  Fourier  transform  of  J(x,y),  which  is  defined  as 


fia,p) 


oo 

-LL 


f(jc,y)  e 


dx  dy 


(2.6) 


from  Eq,  (2.4), 


9  ” 

j(a,)3 )  =  jo(a)  ^  5{a  -  aj  5()3  -  ky), 

m  =  -oo 

(2.7a) 

27an 

(2.7b) 

~  ^  ■> 

where  Jo(tt),  with  a  single  argument,  denotes  a  Fourier  transform  with  respect  to  x  only, 
and  the  tilde  is  used  to  denote  a  Fourier-transformed  quantity.  Equation  (2.7a) 


illustrates  a  common  feature  of  all  periodic  surfaces:  when  a  periodic  surface  is  illuminated 
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by  an  incident  plane  wave,  the  spectrum  of  its  induced  current  is  discrete.  As  a 
consequence,  the  only  scattered  fields  above  and  hplow  the  surface  that  satisfy  the 
boundary  conditions  in  the  plane  of  the  strips  are  plane  waves  whose  transverse  wave 
numbers  are  equal  to  and  ky.  Note  that  we  have  incidentally  derived  the  Floquet  space 
harmonics  -  an  orthogonal  set  of  plane-wave  fields  whose  transverse  phase  constants 
are  given  by  Eq.  (2.7b). 

2,2  Finite  Strip  Grating 

Much  of  what  has  already  been  said  about  the  periodic  strip  grating  applies  equally 
well  to  the  finite  strip  grating.  The  emphasis,  therefore,  will  be  on  the  main  differences 
between  the  two  problems.  The  geometry  is  illustrated  in  Fig.  2.1,  where  a  finite  number 
of  strips  is  assumed  to  exist  at  regular  intervals  along  the  x  axis.  In  this  case,  the 
representation  describing  the  relationship  between  J(x,y)  and  is  trivial,  but  it  is  still 

necessary  to  observe  the  distinction  between  the  two  functions.  For  the  plane  wave 
described  by  Eqs.  (2.1)  and  (2.2), 

J(x,y)  =  JqCx)  ,  *  (2.8) 

where  Jo(Jt)  is  the  current  density  on  the  entire  surface  along  the  cut  y  =  0.  Unlike  the 
previous  case,  the  currents  flowing  on  the  strips  are  no  longer  related  by  a  simple  phase 
shift,  and  in  general,  the  current  on  each  strip  must  be  treated  as  unknown. 

Taking  the  Fourier  transform  of  Eq.  (2.8)  leads  to 

J(a,/3 )  =  2;r  Jo(a)  5(/3  - /ty),  (2.9) 

where,  as  before,  Jo(a)  with  a  single  argument  denotes  a  Fourier  transform  with  respect  to 
X  only. 
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2-3  Periodic  Two-Dimensional  Structures 

Tn  Section  2.1,  a  one-dimensional  structure  was  defined  to  be  one  in  which  the 
surface  geometry  is  independent  of  one  of  the  surface  coordinates.  Next,  the  treatment  is 
extended  to  include  two-dimensional  surfaces,  i.e.,  those,  such  as  the  one  shown  in  Fig. 
2.2,  which  are  functions  of  both  of  the  surface  coordinates. 


y  =  (t  -  m  Ttj  j)  tan  Q 

Figure  2.2.  A  two-dimensional  periodic  surface. 

In  the  general  case,  the  surface  is  considered  to  be  periodic  in  skewed  coordinates, 
with  periodicities  and  Trj2  along  the  ri\  and  772  axes,  respectively.  A  unit  cell  can  be 
defined  as  the  region  -Tx/2<x<  TJ2,  -Tyl2  ^y  ^  Tyl2,  as  in  the  simple  brick  arrangement 
of  Fig.  2.3.  The  unit  cell  is  the  minimum  area  of  the  FSS  that  can  be  reproduced  in  a 
periodic  fashion  to  generate  the  entire  FSS,  including  both  apenure  and  patch  regions 
inside  the  rectangular  boundary. 

The  unit  cell  is  not  unique  and  may  be  chosen  to  suit  the  particular  geometry  of 
interest.  However,  it  is  not  always  possible  to  choose  the  unit  cell  such  that  exactly  one 
element  is  enclosed  within  its  boundaries.  The  two  surfaces  shown  in  Fig.  2.4  are 
common  examples  of  this  situation;  parts  of  at  least  three  elements  must  be  enclosed  within 
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the  unit  cell,  and  care  must  be  taken  to  ensure  that  the  currents  are  continuous  across  these 
fictitious  mathematical  boundaries. 


Figure  2.3.  Simple  brick  arrangement  where  each  unit  cell  contains  exactly 
one  element. 


Figure  2.4.  Common  geometries  for  which  the  unit  cell  encloses  parts  of  at 
least  three  elements. 


Although  the  concept  of  a  unit  cell  is  general  enough  to  describe  any  periodic  surface 
of  the  type  shown  in  Fig.  2.2,  it  is  much  too  restrictive  a  notion  in  practice.  There  is 
nothing  inherent  to  the  problem  that  requires  defining  a  unit  cell.  For  isolated  elements, 
such  as  those  shown  in  Figs.  2.2  -  2.4,  it  is  always  possible  to  represent  the  current 
density  on  the  entire  surface  in  terms  of  the  current  on  a  single  element  of  the  array.  In  this 
case,  the  terminology  “unit  cell”  is  sometimes  used  (very  loosely)  to  denote  a  single 
element.  The  context  will  usually  indicate  the  meaning. 
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It  is  always  possible  to  choose  the  x  axis  of  a  Cartesian  coordinate  system  along  the 
rii  axis  of  the  FSS.  Then  the  diagonal  lines  in  Fig.  2.2,  m  =  0,  ±1,  ±2,  etc.,  are  given  by 

y  =  (jc  -  niT-q^ )  tan  . 

Assuming  that  the  incident  field  is  a  plane  wave  described  by  Eqs.  (2.1)-(2.2),  the  induced 
current  density  J(x,y)  on  the  entire  surface  can  be  written  in  the  form. 


JU,y)  =  JoCx.y)  * 


oo  oo 

^  Six  -y  cotQ  -  ^  Siy  -  nT^^sinQ  ■ , 

m  -  n  = 


(2.10) 


where  JofJt.y)  is  the  current  density  on  a  single  element  of  the  surface.  Except  for  the 
phase-shift  terms,  the  quantity  in  braces  is  a  two-dimensional  array  of  Dirac  delta  functions 
located  at  the  points  of  intersection  of  the  lines  m,n  =  0,  ±1,  ±2,  .  .  .  >»  in  Fig.  2.2. 
Hence,  Eq.  (2. 10)  indicates  that  the  current  density  on  the  entire  surface  can  be  constructed 
by  convolving  the  fundamental  unknown  Jo(x,y)  with  a  two-dimensional  array  of  delta 
functions,  appropriately  modified  by  the  phase  shift  of  the  incident  field. 

As  shown  in  Section  2.1,  this  type  of  representation  is  not  unique.  The  phase  shift  of 
the  incident  field  could  have  been  factored  out  beforehand,  as  follows: 


3ix,y)  = 


oo  oo 

^  Six  -y  cotQ  -  rnTy^)  ^  Siy  -  nT 

m  -  -oo  n  =  -oo 


(2.11) 


Because  of  the  difference  between  the  two  representations,  (2.10)  and  (2.11),  the  two 
currents  Jo(^,y)  are  different  as  well.  Under  certain  conditions,  there  may  be  advantages  to 
using  the  second  type  of  representation,  viz.,  for  electrically  long  elements  or  for  currents 
that  flow  outside  of  unit-cell  boundaries.  However,  the  first  type  of  representation  is  more 
convenient  from  the  point  of  view  of  reconstructing  J(x,y)  from  Jo(j:,y).  The  current  on 
the  patch  at  the  origin  is  equal  to  Jo(J^.j),  and  the  current  on  any  other  patch  is  equal  to 
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JoCr,}')  times  a  constant.  Neither  of  these  statements  is  true  for  the  second  type  of 
representation,  as  is  evident  from  Eq.  (2.11). 

Anticipating  that  the  Fourier  transform  of  J(;c,y)  will  be  needed  later,  the  Fourier 
transform  is  defined  as 

oo  oo 

f(a,)3  )  = /*  f  f{x,y)  dx  dy  .  (2.12) 

«/  — oo  J  — oo 

The  Fourier  transform  of  Eq.  (2. 10)  is  therefore 


J  (a,)3)  =  Jo(a,/3) 


T„  T„  sin^2 

Vi  V2  , 


S  S  S(a-aJ 


(2.13) 


m  =  -oo  n  =  -oo 


where 


(2.14a) 

(2.14b) 


Again,  the  spectrum  of  the  current  is  discrete  and  is  represented  mathematically  as  a  two- 
dimensional  sequence  of  delta  functions  whose  locations  are  the  phase  constants  and 
Pmn  of  the  Floquet  space  harmonics.  The  only  scattered  fields  that  can  exist  above  and 


below  the  surface  are  plane  waves  whose  transverse  wave  numbers  are  equal  to 
am  and  pmn- 


2.4  Finite  Two-Dimensional  Surfaces 

In  progressing  from  periodic  to  finite  frequency  selective  surfaces,  the  problem 
becomes  much  more  difficult  computationally.  For  a  periodic  surface,  the  total  current 
density  on  the  surface  can  be  represented  in  terms  of  the  current  density  on  a  single  unit 
cell;  thus,  it  is  only  necessary  to  work  with  the  equivalent  of  a  single  element  of  the  array. 
In  addition,  since  the  spectrum  of  the  current  is  discrete,  the  moment-method  matrix 
generated  by  applying  the  spectral-Galerkin  technique  is  filled  by  summing  samples  of  the 
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Green's  function,  basis,  and  testing  functions,  in  the  Fourier-transform  domain,  at 
locations  corresponding  to  the  transverse  wave  numbers  of  the  Floquet  space  harmonics. 

In  contrast  to  this,  for  a  surface  which  is  finite  in  two  dimensions,  the  currents  on 
each  cell  are,  in  general,  different,  and  one  is  forced  to  treat  the  current  on  each  individual 
cell  as  unknown.  In  the  context  of  the  method  of  moments,  the  number  of  unknown 
weight  coefficients  for  an  N  X  /V-element  FSS  is  times  the  number  of  coefficients  in  the 
periodic  problem  for  the  same  accuracy.  In  addition  to  being  more  numerous,  the  elements 
of  the  matrix  are  also  more  difficult  to  compute.  The  spectrum  of  the  current  is  now 
continuous  rather  than  discrete,  and  the  moment  matrix  must  be  computed  by  numerically 
integrating  the  product  of  the  Green's  function,  and  basis  and  testing  functions,  in  the 
spectral  domain,  over  the  infinite  a-p  plane. 

As  an  intermediate  step  between  the  periodic  and  the  two-dimensional  finite  surface,  it 
is  expedient  to  consider  first  a  surface  which  is  finite  only  in  one  dimension,  permitting  the 
surface  to  remain  periodic  in  the  other.  In  this  case,  a  finite  number  of  uniformly-spaced 
elements  along  the  jc  axis  are  replicated  in  a  periodic  fashion  along  the  line  y  =  x  tanf2  (see 
Fig.  2.2)*  The  currents  on  a  given  row  of  elements  are  no  longer  simply  related;  however, 
assuming  that  the  incident  field  is  a  plane  wave,  the  currents  on  any  row  of  elements  are 
related  to  the  currents  on  any  other  row  through  the  phase  shift  of  the  incident  field.  As  a 
result,  if  the  incident  field  is  derived  from  the  electric  and  magnetic  vector  potentials  (2. 1  )- 
(2.2),  then  the  induced  current  density  on  the  entire  surface  can  be  represented  by  a 
convolution  of  Jo(x,y),  the  current  density  on  a  unit  row  of  elements,  and  an  infinite 
periodic  sequence  of  Dirac  delta  functions,  appropriately  phase-shifted  to  account  for  the 
effect  of  the  incident  field: 


JCc,y)  =  Jo(x,y)  ^  S(y  -  sinQ )  Six  -ycot{2 )  e'^y^ 


.  (2.15) 
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The  sequence  of  delta  functions  used  here  to  construct  the  total  current  from  JoU.y)  is  the 
sequence  of  delta  functions  in  Fig.  2.2  corresponding  to  m  =  0. 

Taking  the  Fourier  transform  of  (2.15)  leads  directly  to 


j(a.j3)  =  jo(tt,/?)  ^  Sip-j5,{a)),  (2.16) 

sin^2  n  =  — oo 

where 

Except  for  the  special  case  =  90°,  the  phase  constants  Pn  will  be  functions  of  a. 

As  stated  previously,  for  a  frequency  selective  surface  which  is  finite  in  two 
dimensions,  there  is  no  general  relationship  between  the  current  on  the  entire  surface  and 
the  current  on  any  pan  of  the  surface.  In  the  notation  of  this  section, 

J(x,y)  =  JoU,y), 

where  Jo(.^.y)  is  the  current  density  on  the  entire  two-dimensional  surface. 

In  summary,  the  Fourier  transform  of  the  total  induced  current  density  J(a,/3)  can  be 
expressed  as  the  product  of  Joicc,p),  the  Fourier  transform  of  the  current  density  on  a 
subsection  of  the  surface,  times  a  periodicity  function  ri(a,p),  which  is  either  a  constant, 
or  a  one-  or  two-dimensional  array  of  Dirac  delta  functions. 


n{a,p) . 


(2.18) 
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CHAPTER  3 

THE  SPECTRAL-GALERKIN  TECHNIQUE 


3ol  Derivation  of  the  Spectral-Domain  Green's  Function 

One  of  the  main  reasons  for  choosing  the  spectral-domain  formulation  is  the  ability  to 
derive  the  spectral-domain  Green’s  function  -  a  transform  relationship  between  the 
scattered  electric  field  in  any  transverse  plane,  due  to  the  surface  current  J  in  any  other 
plane. 


G^y(a,p)  G^y(a,p) 


Ey’(a,fi) 


Gy^{a,fi)  Gyy(a,p) 


Jy<a,p) 


(3.1) 


Here  we  have  followed  the  convention  of  using  the  tilde  to  indicate  a  Fourier-transformed 
quantity,  with  the  Fourier  transform  defined  as 


f{a,p  )=j  j  i{x,y  )  e  e  dx  dy  . 

Interpreting  the  elements  of  G  from  Eq.  (3.1),  Gx'x,  for  example,  is  the  Fourier  transform 
of  Ex',  the  x'  component  of  the  scattered  electric  field  in  the  plane  z'  =  0,  divided  by  the 
Fourier  transform  of  an  infinite  planar  current  sheet  Jxix^y")  =  e/(o*’  Py'). 

There  are  very  few  practical  situations  in  which  the  printed  surface  can  be  considered 
as  free-standing  -  in  the  absence  of  any  supporting  dielectric  layers.  Most  often,  the 
properties  of  the  layers  are  fixed  by  mechanical  constraints,  i.e.,  because  of  weight 
limitations,  or  requirements  for  temperature  stability,  shape  and  rigidity,  FSS  support 
structures  constructed  with  alternate  layers  of  honeycomb  and  thin  film  are  common  in  FSS 
applications.  Optimization  techniques  have  also  been  applied  in  order  to  design  dielectric 
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layers  to  modify  the  electromagnetic  properties  of  frequency  selective  surfaces  [14],  The 
interest  in  free-standing  surfaces  is  therefore  limited,  and  it  is  essential  that  the  FSS  be 
considered  in  an  environment  consisting  of  a  number  of  arbitrary  dielectric  layers. 

The  usual  method  of  incorporating  the  effects  of  the  dielectric  support  is  the 
scattering-matrix  approach  [4,  15,  16].  If  a  plane-wave  field  is  incident  on  a  periodic 
surface,  the  scattered  field  is  composed  of  an  infinite  number  of  discrete  plane  waves, 
whose  wave  numbers  am,  Pmn  the  transverse  phase  constants  of  the  Floquet  space 
harmonics.  If  any  one  of  these  plane-wave  fields  is  now  incident  upon  another  surface 
having  the  same  periodicity,  the  scattered  fields  will  consist  of  the  same  set  of  Floquet 
harmonics.  In  general,  then,  if  a  multilayered  structure  is  composed  only  of  uniform 
homogeneous  dielectric  layers  and  printed  periodic  surfaces  having  the  same  lattice  {T-qi, 
T-q2,  •^2)>  then  the  fields  existing  throughout  the  region  can  be  represented  exactly  by  a 
doubly-infinite  set  of  plane-wave  fields  called  the  the  Floquet  space  harmonics.  In  a 
lossless  medium,  a  finite  number  of  these  plane-wave  fields  will  be  propagating,  and  the 
plane-wave  fields  corresponding  to  higher  values  of  m  and  n  will  eventually  be 
insignificant  in  comparison  to  the  dominant  terms.  It  is  therefore  only  necessary  to 
consider  a  finite  number  of  Floquet  harmonics  in  order  to  obtain  an  accurate  representation 
for  the  fields  in  each  layer. 

Scattering  from  a  given  layer  can  now  be  described  in  terms  of  a  matrix,  whose 
elements  relate  the  amount  of  power  in  a  given  scattered  plane-wave  field  due  to  an  incident 
plane  wave  from  another  direction.  A  scattering  matrix  is  computed  individually  for  each 
layer,  and  the  matrices  are  cascaded  to  predict  the  scattered  fields  from  the  composite 
structure. 

Because  there  is  nothing  analogous  to  the  scattenng-matrix  approach  for  cascading 
finite  frequency  selective  surfaces,  it  is  desirable  to  develop  the  spectral-domain  Green's 
function  from  a  point  of  view  which  is  general  enough  to  be  able  to  handle  the  geometry 
shown  in  Fig.  3.1,  which  is  typical  in  practice  [14,  17].  Here,  it  is  assumed  that  a  single 
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layer  of  metallization  is  embedded  in  an  /V-layer  stratification,  where  each  material  is 
assumed  to  be  an  infinite,  uniform,  homogeneous  dielectric  layer  of  complex  relative 
dielectnc  permittivity  £r  and  magnetic  permeability  fj^. 
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Figure  3.1.  Typical  dielectric  support  used  in  FSS  applications. 


The  derivation  of  the  spectral  Green's  function  is  substantially  the  method  of  Itoh 
[18],  which  is  based  on  the  transmission-line  analogy  for  a  plane  stratified  medium.  For 
completeness,  and  for  the  benefit  of  a  different  point  of  view,  a  detailed  derivation  of  the 
spectral  Green's  function  is  provided  below. 

The  transmission-line  analogy  is  based  on  the  following  field  properties;  (1)  any 
arbitrary  field  in  a  homogeneous  source-free  region  can  be  expressed  as  the  sum  of  a  TE 
and  a  TM  field;  and  (2)  a  plane  wave  propagating  in  a  plane  stratified  medium,  polarized 
either  TE  or  TM  with  respect  to  the  surface  normal,  will  remain  TE  or  TM  throughout  each 
layer  of  the  stratification.  Because  of  the  similarity  between  the  two-dimensional  Fourier 
transforms  of  Maxwell's  equations  and  the  transmission-line  equations,  for  each 
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polarization,  a  length  of  transmission  line  can  be  assigned  to  each  layer  of  the  structure  and 
the  transmission-line  theory  can  be  used  to  derive  the  Green's  function. 

In  general,  an  arbitrary  current  sheet  Jx’i.x'.y')  orJy’ix'.y')  =  produces 

both  TE  and  TM  fields.  Therefore,  in  order  to  apply  the  equivalent  transmission-line 
analogy,  the  fields  must  be  separated  into  their  TE  and  TM  components.  To  do  this,  it  is 
advantageous  to  perform  a  rotation  of  the  coordinate  axes  [18]  in  order  to  simplify  the 
association  between  the  field  quantities  E  and  H  and  the  scalar  quantities  V  and  /. 

The  particular  rotation  of  coordinates  that  achieves  this  simplification  is  determined  by 
the  following  considerations:  in  order  for  the  phases  of  the  fields  to  match  at  every  point 
on  the  boundary  planes,  the  fields  set  up  in  each  region  due  to  the  current  sheet  must  be 
plane  waves.  Furthermore,  the  phase  constants  a  and  p  of  the  assumed  current 
distribution  must  be  the  x  and  y  components,  respectively,  of  the  wave  vectors  of  those 
plane-wave  fields.  The  desired  simplification  results  if  the  coordinate  system  is  rotated 
about  the  2'  axis  until  TE  and  TM  with  respect  to  z'  are  also  TE  and  TM  with  respect  to  the 
plane  of  the  wave  vectors.  In  the  new  system,  both  E  and  H  will  have  exactly  one 
tangential  component  at  the  boundary  planes,  regardless  of  polarization.  Equally  as 
important,  this  coordinate  transformation  has  the  property  that  when  the  current  sheet  Jx'  or 
Jy'  is  represented  in  the  new  coordinate  system,  one  of  its  components  will  give  rise  to  a 
field  which  is  TE  to  z'  and  the  other  component  will  give  rise  to  a  field  which  is  TM  to  z'. 

The  rotation  of  coordinates  obtained  by  rotating  the  original  through  an  angle 

(p  about  the  z'  axis  satisfies  mU  of  the  conditions  necessary  to  make  a  simple  transition 
between  the  fields  problem  and  the  transmission- line  problem.  In  the  next  section, 
however,  it  will  be  desirable  to  have  the  incident  field  propagating  in  the  +z  direction  to 
agree  with  the  usual  transmission-line  conventions.  Therefore,  the  final  coordinate  system 
is  obtained  by  reversing  the  z'  axis.  The  result  is  shown  in  Fig.  3.2,  where,  as  long  as  it  is 
necessary  to  make  a  distinction  between  the  two,  the  original  coordinate  system  will  be 
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designated  by  the  primed  coordinates  x',  y',  z  and  the  new  coordinate  system  by  the 


unprimed  coordinates  x,  y,  z; 


r sin0  ir  X' 


yj  Lsin<^  -  COS<^JLy 


fcos^  sin^^lfi 


y’J  [sin^ -cosily 


z=  -  z 


Figure  3.2.  Coordinate  transformation  defined  for  the  equivalent 
transmission-line  analogy. 


An  arbitrary  A-layer  medium  is  shown  in  Fig.  3.3,  with  its  corresponding  equivalent 
transmission-line  system,  where  the  polarization  p  refers  to  TE  or  TM.  In  the  layer,  if 
the  z-dependence  of  the  fields  traveling  in  the  +z  direction  is  of  the  form  e'Ti^ ,  then  the 
propagation  constants  yi  are  given  by 

7i  ='J a^+  P  ^  -  Erikl  ,  (3.3) 
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where  is  the  free-space  wave  number.  The  characteristic  impedances  of  the  equivalent 
transmission  lines  are  defined  as  the  wave  impedances  of  the  plane  waves  propagating  (or 
attenuating)  in  the  +z  direction  in  each  medium.  From  Maxwell's  curl  equations, 


7TE  ~^y 

~  H,-  y, 

(3.4a) 

II 

II 

(3.4b) 

With  the  equivalent  transmission  lines  defined  by  (3.3)  and  (3.4),  the  problem  can  now  be 
solved  using  a  handful  of  highly-developed,  systematic  techniques. 


X' 


Figure  3.3.  /V-layer  dielectric  support  and  its  equivalent  transmission-line 
network. 
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In  a  given  length  of  transmission  line,  the  general  solution  to  the  telegrapher's 
equations  can  be  expressed  as  a  superposition  of  a  +z-directed  and  a  -z-directed  wave, 

V(z)  =  V  V~e'^=V  \z)  +  V  "(z) 

/(z) 

where  V+  and  V-  are  arbitrary  constants.  If  the  generalized  reflection  coefficient  Ffz)  is 
defined  as  the  reflected-wave  voltage  divided  by  the  incident-wave  voltage,  i.e., 


r(z)  = 


v~iz) 

v^iz) 


then  the  voltage  and  current  on  the  line  can  be  expressed  in  the  form. 


V(z)  =  I^(z)[l -t-Hz)]  (3.5a) 

/(z)  =  ^V^(z)[l-r(z)]  .  (3.5b) 


Taking  the  ratio  of  the  voltage  and  current  leads  to  an  expression  for  the  impedance  Z(z), 


Viz)  ^  l-^Hz) 
/(z)  "^^l-Hz)  ’ 


(3.6) 


or,  rewriting  this  in  terms  of  r(z), 


r(z)  = 


Z.(2)-l 

Z,(z)  +  1 


where  Z;,(z)  is  the  normalized  impedance  Z(z)/Zo. 
it  follows  that 


(3.7) 

Since  V+(z)  =  V'^e  l^  and  V '(z  )  = 


r(z,)  =  r(z2) 


(3,8) 


for  any  two  points  z\  and  Z2  on  a  section  of  uniform  transmission  line. 

This  brief  review  suggests  the  following  strategy  [19];  working  from  the  load  end  to 
the  source  end  of  the  transmission  line,  in  any  uniform  section  of  line,  (1)  normalize  the 
impedance  at  the  load  end  by  the  characteristic  impedance  of  the  section,  (2)  calculate  the 
reflection  coefficient  at  the  load  end,  (3)  rotate  the  reflection  coefficient  toward  the  source 
using  Eq.  (3.8),  (4)  compute  the  normalized  impedance  at  the  source  end,  (5)  unnormalize, 
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and  (6)  advance  to  the  next  section.  The  operations  are  highly  repetitive  and  are  efficiently 
implemented  on  a  digital  computer. 

Once  the  input  impedance  seen  by  the  source  has  been  determined  by  the  above 
procedure,  the  remaining  steps  are  independent  of  the  specifics  of  the  problem.  The 
boundary  conditions  satisfied  by  the  fields  E(z)  and  H(z)  at  z  =  0  are 


zxlE(0l-E(0')]  =  0  (3.9a) 

zx[H(0-")-H(0")]=J,  (3.9b) 

or,  expanding  Eq.  (3.9)  in  Cartesian  coordinates, 

=  (3.10a) 

/;,  =  -(  HyiO^)  -  Hy{0~) ) .  (3. 10b) 


Assuming  for  the  moment  J(x',y')  =  + /^y’),  then,  according  to  the  coordinate 

transformation, 

J  =  [  xcos^  +  y  sin^]  (3.11) 

The  boundary  condition  satisfied  by  the  magnetic  field  at  the  source  plane  is  that  the 
discontinuity  in  Htan  is  equal  to  the  surface  current  density  J.  Since  the  tangential  magnetic 
field  is  analogous  to  the  current,  symbolized  by  the  notation  Htan  <->  /,  a  discontinuity  in 
Htan  corresponds  to  a  shunt  current  source  in  the  equivalent  transmission  line  model.  For 
the  TE  fields  (Ey,  Hx,  Hz)  and  their  corresponding  transmission-line  system,  the 
appropriate  boundary  condition  on  Htan  and  Kirchoff  s  current  law  at  the  source  are  (Fig. 
3.4) 

KCL:  /,  =  /(0'")-[-/(0-)]  (3.12) 

Jy  =  HjSi")-HjSy) .  (3.13) 

Comparing  (3.11),  (3.12),  and  (3.13),  the  value  of  the  equivalent  current  source  I  s=  sintp, 
and,  since  V  =  IsZ^^  and  -V  (Eq.  (3.4)), 

Ef  =  sin0  =  e]1 


(3.14) 
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I  (O'  I 


I  (0^) 


'.0 


:z,^  z^y/Zp 


Figure  3.4.  Equivalent  circuits  at  the  source  terminals  for  the  transmission- 
line  system  of  Fig.  3.3. 


For  TM  polarization  (Hy,  Ex,  Ez), 


KCL:  4  =  /(0'')-[-/(0  )] 
BC:  -4  =  //y(0'")-//y(0') 


(3.15) 

(3.16) 


Again,  comparing  (3.1 1),  (3.15),  and  (3.16),  Is  -  -cos(p,  and,  since  V  = 


^  +F(Eq.  (3.4)), 


£™  =  -z™  COS.? 


(3.17) 


In  the  remaining  steps,  the  fields  are  converted  back  to  the  original  coordinate  system,  and 
the  TE  and  TM  components  are  superposed  to  obtain  the  total  field. 


£,.  =  [  £"sin?  +  £™  cos«  1  c'»'  =  -  ( Z^f  sinV  +  Z™  cosV  1  c'”'  (3.18a) 

£,.  =  I  -£^cos<.  +  eIS*  sin?  1  =  +  ( z[„^  -  Z™  1  sin?  cos?  c““'  (3. 18b) 


(3.18a) 


or,  from  the  definition  of  G  (Eq.  (3. 1 )), 


_  r  "tTE  •  2  .  .  .yTM  1 

G,'x'  =  —  =  -  I  Zi„  Sin  0  +  Z,„  cos  0  J 


(3.19a) 


Gv';,'  =  —  =  +  [  ]  sin0  COS0 


(3.19b) 
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Sincp  and  cos(p  are  related  to  a  and  /3  by  the  following  relations: 


sin^= - 

(3.20a) 

costp  - - 

(3.20b) 

as  is  evident  from  Fig.  3.2.  A  similar  procedure  can  be 

Gx'yicc,p )  and  Gyyia,^ ).  Only  the  results  are  summarized  below: 

used  to  derive 

G^y  =  ^  =  -H  [  Zyf  -  Z™  ]  sm<p  cos(p 

jy 

(3.21a) 

Gyy  =  ^  =  -  [  z'jj  COS^0  -H  Z™  sinV  ]  • 

(3.21b) 

As  an  important  special  case  of  the  above,  the  free-space  Green's  function  can  be  written 
immediately  by  recognizing  that 


zTf=z^//zI^ 


jm 

2/0 


yTM 

^in 


_  ^TM  ,,  ^TM 

-■^00  //•^oo 


To 

Ijae^  ’ 


(3.22a) 

(3.22b) 


TE  TM 

where  and  Zqq  are  the  characteristic  impedances  for  the  TE  and  TM  plane-wave  fields 
propagating  in  free-space,  and  Yq  is  the  free-space  propagation  constant,  as  defined  by  Eq. 
(3.3)  with  =  1.  From  (3.19)  and  (3.21), 


G{a,P)  = 


1 

2/(y€b7o 


^0  —  —ccp 

-ap  li-f 


(3.23) 


As  a  second  example,  the  Green's  function  for  a  two-layer  dielectric  suppon  is 
derived,  where  the  current  sheet  is  assumed  to  be  embedded  between  two  layers 
characterized  by  complex  relative  dielectric  permittivities  Erx,  £n  and  thicknesses  t\,  t2  as 
shown  in  Fig.  3.5.  When  the  two  transmission  lines  extending  to  z  =  ±<»  are  replaced  by 
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their  characteristic  impedances,  the  impedances  and  as  viewed  in  the  +z  and  -z 
directions,  respectively,  are  written. 


^02  ^00^0^^72/2 

Z02 - 

^00  +  2o20oth72r2 


^01  +  ^ooooth7iti 
^01  • 
^00  ^oi0oth7iri 


The  total  input  impedance  seen  by  the  source  is  the  parallel  combination  of  Z^  and  Z_, 

m  y  y 

ZoiZ^2[2o2  +  2a)C0th72r2  ][Zoi  +  Z^cothyiti  ] 

^02^2^00  ■•■  ^oioothyitj  ][Zo2  +  Z^coth72r2  ]  +  ZqiIZqi  +  Z^cothyitj  ]{Zqq  +  Zo2coth72i2  ] 

(3.24) 

The  Green's  function  is  then  given  by  Eqs.  (3.19)  and  (3.21). 


1  + 


Zr  _  'yt' 

00 


z* 


=  zs 


Znn  +  Zn 


£1^-272^2 


02‘ 


1  - 


ZP  _  7P 
00 


ZS0  +  Z02 


£2^-2y,r2 


1  + 


ZP  ^  7P 
00 


2;  =  zg, 


2oO  +  2gi 


2ig-2r,h 


1  - 


ZP  _  yP 
00 


4o  +  4i 


£1^-27, 


3.2  Asymptotic  Behavior  of  the  Green's  Function 


Asymptotic  expressions  for  the  spectral  Green's  function  are  easily  found.  If  » 
\£ri\k^  and  ^pT,  >  n  for  each  finite-thickness  layer,  then  asymptotically 

Yi  kp  (3.25) 

coth7/i— >1.  (3.26) 


The  characteristic  impedances  of  the  equivalent  transmission  lines  are,  therefore, 


.7TE 

^Oi 


jjm 

k„ 


(3.27a) 


^TM 

^Oi 


jOJCriEo 


(3.27b) 
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Figure  3.5.  Geometry  of  two-layer  dielectric  support  and  its  equivalent 
transmission-line  system. 


According  to  Eq.  (3.27),  the  characteristic  impedances  for  the  TE  transmission-line  system 
are  independent  of  the  parameters  of  the  layers.  The  characteristic  impedances  of  the  TM 
transmission  line  system  differ  only  by  a  multiplicative  constant.  As  a  result,  the  input 
impedance  zj^  is  the  parallel  combination  of  two  terms  of  the  form  jcofjJkp, 


(3.28) 


Furthermore,  it  is  not  difficult  to  verify  that  under  the  approximations  (3.24)-(3.27)  the 
TM 

input  impedance  Z-^  depends  only  on  the  dielectric  constants  of  the  layers  nearest  the 


source. 
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jTM 

'  "  jQ}£o(£^  +  £,J  ' 

where  and  fr-  the  relative  permittivities  of  the  two  layers  adjacent  to 
sheet.  Substituting  the  asymptotic  formulas  (3.28)  and  (3.29)  into  Eqs.  (3.19) 
the  Green's  function  becomes 


(3.29) 

the  current 
and  (3.21), 


-1 

.  a' 

.  2(a^  +  /3^)  (£r+  +  £r-)_ 

-ap 

kl  1  ' 

-1 

-kla^ 

.  2(a^  +  j3^)  (^+  +  er-)_ 

(3.30a) 

(3.30b) 

(3.30c) 


These  results  are  useful  for  several  reasons.  First,  since  the  Green’s  function  appears 
in  the  product  j,t(o',/J)G(a,/?)jy(a,/3),  the  asymptotic  form  of  the  Green's  function 

determines  the  minimum  rate  of  decay  necessary  in  the  Fourier  transforms  of  the  basis  and 
testing  functions  in  order  for  the  inner  products  to  converge.  If,  in  polar  coordinates,  the 
product  does  not  decay  faster  than  l/kp,  the  inner  product  will  not  converge,  and  a 

smoother  set  of  basis  and  testing  functions  must  be  chosen.  Second,  since  the  asymptotic 
forms  are  considerably  simpler  than  the  general  forms,  it  is  possible,  for  a  judicious  choice 
of  basis  and  testing  functions,  to  perform  the  asymptotic  integrations  exactly  [20].  This 
would  offer  advantages  in  terms  of  both  speed  and  accuracy.  However,  even  if  the  inner 
products  cannot  be  evaluated  exactly  using  the  asymptotic  form  of  the  Green’s  function,  the 
asymptotic  form  requires  considerably  less  computational  effort  than  the  general  form. 
Numerical  experiments  have  shown  that  significant  savings  in  computation  time  can  be 
achieved  by  using  the  asymptotic  form  when  kp  becomes  sufficiently  large. 


3.3  Singularities  of  the  Spectral  Green's  Function 


The  spectral-domain  Green's  function  for  sources  in  a  lossless  plane  stratified 

medium  has  singularities  of  two  different  types.  First,  since  the  propagation  constants  yi 

are  defined  in  terms  of  the  square  root  function,  which  is  a  multivalued  function,  there  are 

branch  points  at  the  locations  where  kp'^  =  these  points  the 

derivative  of  the  Green's  function  is  undefined.  There  are  also  nonintegrable  singularities 

caused  by  zeros  in  the  denominators  of  the  input  impedances  and  Z^^  which 

correspond  to  guided  modes  in  the  dielectric-slab  waveguide. 

Beyond  this,  very  little  can  be  said  about  the  singularities  without  considering  a 

specific  example.  Returning  to  the  two- layer  example  discussed  in  Sec.  3.1,  if  the 

substrate  and  superstrate  are  identical,  such  that  £ri  =  Erl  =  =  ^2  =  ? »  and  yi  =  yi  =  y, 

TE  TIvl 

then  the  denominators  Aje,  zixM  of  the  impedances  Z  -^  and  Z  -^  ,  respectively,  are  given 
by 

(-JcoEof 

Aje  =  — T7~~  '^^ry(^rYo  +  }coth7r)(e^7ocoth7r  +y)  (3.31a) 

W 

^TM  =  2)<r+  rocoth7r)()coth7r  4  yg) .  (3.31b) 

If  pole-zero  cancellation  does  not  occur,  then  the  non-branch-type  singularities  of  the 


Green's  function  will  be  solutions  of 

(^rTo  +  Tcothyr)  =  0  (3.32a) 

(e,7ocoth7r  +  y)  =  0  (3.32b) 

{y-¥  7ocoth7r)  =  0  (3.32c) 

(}Coth7r  -I-  7o)  =  0  .  (3.32d) 


Since  Tj  is  a  function  of  k^^  only,  i.e.,  since  yi  is  not  a  function  of  the  polar  angle,  the 
singularities  of  Zj^  andZ^^  lie  on  concentric  circles  about  the  origin  in  the  a-P  plane. 
Since  Eqs.  f3.32a-d)  are  transcendental  equations,  the  locations  of  the  singularities  can 
only  be  determined  numerically.  However,  if  the  dielectrics  are  lossless,  the  following 
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general  observations  can  be  made.  First,  the  number  and  locations  of  the  singularities 
depend  on  frequency,  the  relative  dielectric  constants,  and  thicknesses  of  the  layers.  If  any 
two  of  these  parameters  are  held  constant  while  the  other  is  varied,  higher  frequencies, 
higher  dielectric  constants  and  greater  thicknesses  all  tend  to  produce  more  singularities  in 


kp  >  where  is  the  maximum  of  ji  is 

real  and  positive,  cothyjr  is  real  and  positive,  and  each  term  in  Eqs.  (3.32a-d)  is  the  sum  of 
two  nonzero  positive  numbers.  In  such  a  case,  Eqs.  (3.32a-d)  have  no  real  solutions. 
Recalling  that  the  branch-type  singularities  occur  at  kp  =  ^£r^ko  ,  there  can  be  no 


singularities  of  either  type  in  the  region  kp  >  '^£rmax^^-  Likewise,  for  kp  <  ko,  a  similar 
argument  holds.  Therefore,  all  of  the  singularities  lie  somewhere  in  the  annular  region  ko  < 
kp  <  Lis^re  3.6  is  a  plot  of  the  Green’s  function  for  a  structure  with  two 

identical  layers,  characterized  by  Er  ~  4.0,  t  =  Ao/6  along  the  radial  line  a  =  /?.  Eventually, 
Gxxicc>P )  and  Gyyia.P )  increase  linearly  with  kp^  becoming  unbounded  at  oo. 


Krho(radlm) 


Gxx 

Gxy 

Gyy 


Figure  3,6.  Spectral  Green's  function  for  the  fields  at  the  center  of  two 
identical  layers  of  permittivity  =  4.0  and  thickness 
t  =  Ao/6  (frequency  =  10  GHz). 
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3.4  The  Spectral-Galerkin  Equation 

As  was  shown  in  Sec.  3.1,  the  general  relationship  between  the  scattered  electric  field 
and  an  arbitrary  surface  current  J  radiating  in  the  presence  of  a  number  of  infinite, 
homogeneous  dielectric  layers  can  be  written  in  the  form 


E^ia,p) 


^Eyia.p) 


G,./a.p)  G,yia,p) 


Gyy{a,p) 


yy 


Jyia.p) 


^Jy{a,p')  _ 


The  current  density  J  is  now  interpreted  as  the  induced  current  on  a  printed  frequency 
selective  surface  (whether  finite  or  periodic),  given  by  one  of  the  representations  (2.4), 
(2.8),  (2.10),  or  (2.15).  Except  for  one  or  two  special  cases  [10],  the  current  Jo(.r,y) 
cannot  be  determined  analytically,  and  numerical  techniques  must  be  used  in  order  to 
compute  it. 

The  spectral-domain  method  of  moments  proceeds  as  follows.  Choosing  a  suitable 
set  of  linearly-independent  basis  functions  {jy},  the  current  Jo  is  expanded  in  the  form, 


•4o 

-•^>0 


Jxj 

Jyj  . 


(3.33) 


where  [cj]  is  a  set  of  complex-valued  weight  coefficients  which  are  yet  to  be  determined. 
There  are  four  main  considerations  involved  in  choosing  an  appropriate  set  of  basis 
functions  in  the  spectral-domain  method  of  moments.  First,  as  with  any  moment-method 
solution,  the  basis  functions  must  be  linearly  independent.  If  any  one  basis  function  could 
be  represented  in  terms  of  a  linear  combination  of  two  other  basis  functions,  for  example, 
then  one  column  of  the  moment  matrix  would  be  representable  in  terms  of  two  others,  and 
the  matrix  would  be  singular.  Second,  since  the  method  of  solution  is  a  transform-domain 
technique,  the  basis  and  testing  functions  must  be  Fourier  transformable.  This  is  not  a 
serious  limitation,  since  all,  or  nearly  all,  basis  functions  used  in  practice  are  Fourier 
transformable.  A  third  consideration,  which  is  especially  significant  for  the  analysis  of 
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finite  frequency  selective  surfaces,  is  that  the  basis  must  be  capable  of  representing  the 
fundamental  unknown  with  a  small  number  of  terms.  This  is  because  the  matrix  size  is 
proportional  to  for  an  /Vx/V-element  array  For  example,  for  the  entire-domain  basis  set 
of  Orta,  Tascone,  and  Zich  [21],  about  twelve  basis  functions  per  arm  are  necessary  to 
accurately  represent  the  current  on  a  particular  tripole  element.  At  36  basis  functions  per 
element,  even  a  small  10x10  array  would  require  filling  and  storing  a  3600x3600-element 
matrix.  A  fourth  consideration  that  must  be  observed  in  choosing  a  basis  is  the  requirement 
that  the  basis  and  testing  functions  together  have  sufficient  smoothness  to  guarantee  that  the 
integrals  and  summations  generated  by  the  spectral  method  of  moments  converge.  It  is  also 
desirable  that  the  basis  functions  model,  at  least  approximately,  the  known  physical 
behavior  of  the  current.  For  example,  the  cunent  flowing  parallel  to  the  edge  of  a  printed 
element  exhibits  an  edge  singularity  proportional  to  where  (x  -  Xq)  is  the 

distance  from  the  edge. 

Substituting  for  J  and  Jo,  the  scattered  electric  field  is  of  the  form 

Next,  an  inner  product  is  defined  satisfying  properties  (l-2)-(l-4)  of  Harrington  [22], 


.  OO  _ 


^xx  ^x  y  Jxj 
Gyj^  Gyy  _  _Jyj. 


ha,^)  .  (3.34) 


(  f,  g  )  =  /  f^(A:)  g{x,y)dx  ,  (stnp  grating)  , 

oo  oo 

(  f'  g  )  =  /  /  Cx,y)  g{x,y)dxdy  ,  (two-dimensional  element) 

oo  •/— oo 


(3.35a) 

(3.35b) 


where  the  superscript  “T’  signifies  transpose.  Choosing  the  testing  functions  to  be  the 
complex  conjugates  of  the  basis  functions  [23],  and  forming  the  inner  product  of  Eq. 
(3.34)  with  each  testing  function  lead  to  a  system  of  /  equations  in  /  unknowns 
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J 

E' )  =  j J  ff 


Gia,p)]j{a,[5)nia,(5}  dadfi , 

(3.36) 


where  the  iagger  signifies  complex  conjugate  transpose.  The  double  integral  over  x 
and  y  can  be  recognized  as  a  Fourier-transform-type  integral,  which  can  be  done 
analytically,  resulting  in  the  equation 
I  OO  OO 


f  f  ji(a,p)  G(a,p)  jj(a,p)n(a,p)  dad(i  E 

Att  J  ^  1 


me 


>,  /■=  1.  2 . / 


(3.37) 


As  a  result  of  testing  with  jj,  the  equation  has  been  restricted  to  the  surfaces  of  the 
conductors.  The  boundary  condition  has  therefore  been  enforced  by  replacing  the  scattered 
field  with  the  negative  of  the  incident  field  on  the  right-hand  side.  The  specific  form  of  Eq. 

(3.37)  will  depend  on  the  specifics  of  the  problem,  but  for  the  five  general  types  of 
problems  discussed  in  Chapter  2,  the  spectral-Galerkin  equation  is  summarized  below. 


1.  Periodic  Strip  Grating  (Fig.  2.1): 


n{a,l3)  =  ^  Y,  dia  -  aj  b{{3  -  ky) ,  + 


(3.38) 


X  Gia^.ky)  j/aJe'V  =  f  jj(x)  E‘"‘'(x,y)  dx  ,  1=1,2,  ...,  I . 

(3.39) 


i  =  \  m- 


2.  Finite  Strip  Grating  (Fig.  2.1): 

rKa,p)  =  27t5{p-ky) ,  (3.40) 

^  f  OO  oo 

J  Gia,ky)  j/a)e'Wo:  =  f  j-(x)  E‘"‘'(x,y)  dx  ,  i  =  1,  2,  ...,  /  . 


(3.41) 


3.  Periodic  Two-Dimensional  Surfaces  (Fig.  2.2): 


oo  oo 


'  ^2  m  =  -oon  =  -°o 

(3.42) 

iTtm 

®/n  ~  T  ’ 

n, 

^  2;m  2;tm  ^ 

“  T„  sinr2  T„  ’ 

"2 

(3.43a) 

(3.43b) 
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Yc 

00  oa 

X  X  ii  G(a;„,^;„„)  j/ 

;n  =  -00  n  =  -00 

oe  00 

=j  j  jj(x,y)  E‘''‘^(j:,y)  dxdy  ,  t  =  1,  2,  ...,  / 

(3.44) 

4,  Finite-Width  Two-Dimensional  Surfaces  (Fig.  2.2): 

OO 

i7(a,/3)  = . . 

T yj  ^sin^2  n  =  —00 

(3.45) 

/5„(a)  =  j  (a  UcotX2  -h  ky , 

^2 

(3.46) 

^  I 

27tr„  sini2^^J 

^2  J=l 

00 

X  /  j?(aA)  G(aA)  j/a-i3„)da 

rt  =  -oo 

00  00 

=  f  f  jl(^,y)  E'’'''(x,y)  dxdy  ,  /  =  1,  2,  ...,  I . 

•^—00  oe 

(3.47) 

5.  Surfaces  That  Are  Finite  in  Two- Dimensions  (Fig.  2.2); 


IXa,P)  =  1  , 


(3.48) 
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T„  Ty.  sin^2 

”2  y  =  1 


ijia.p)  Gia,p)  jj(a,p)dadp 

ee  oo 

=  /*  f  E‘"‘'U,>’)  dxdy  ,  1=1,2, ...,  I .  (3.49) 

tf  00  a/  ae 


From  this  summary,  it  is  possible  to  see  that  the  difficulty  involved  in  computing  the 
moment  matrix  ranges  from  performing  a  summation  over  samples  of  the  Green's  function, 
basis,  and  testing  functions  at  the  discrete  values  a^,  Pmn  to  a  complete  two-dimensional 
integration  over  the  infinte  a-fi  plane. 


3.5  Incident  Field 

The  spectral-Galerkin  equation  derived  in  the  previous  section  is  written  in  terms  of 
the  quantity  defined  as  the  electric  field  reaching  the  plane  of  the  printed  conducting 
surface  with  the  conductors  removed,  i.e.,  in  the  presence  of  the  dielectric  suppon.  For 
plane-wave  incidence,  the  tangential  field  E,"”^  is  of  the  form, 

Er=  (  (3,50) 

where  Ex'o  and  Eyo  are  complex  constants.  When  combined  with  the  testing  functions,  the 
inner  product  (3.28)  becomes 


(j/,  E;'“')  =  JJ  hyix’y)  ] 


EyO 


+  (3.51) 


This  can  be  recognized  as  a  Fourier-transform  type  of  operation,  allowing  the  integration  to 
be  done  analytically,  resulting  in 

(jj,  ED  =  E,.^[  hAkx;ky)  ]*  +  £yo[  hyiK;ky)  ]*  .  (3.52) 

For  the  strip  grating, 

oe 

(jj,  E;"">  =  f  \]{x’)ET‘^{x',y')dx'  =  {£^-0  +  £yo  liy{k,.)]e>’^y^‘ .  (3.53) 

•Y*  =  — OO 


In  order  to  find  the  constants  £x'o  and  Ey^^,  we  will  again  make  use  of  the 
transmission-line  analogy  and  rotation  of  coordinates  defined  in  the  derivation  of  the 
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Green's  function.  In  passing,  we  note  that  the  fields  in  a  plane  stratified  medium  can  also 
be  determined  by  writing  a  representation  for  the  fields  in  each  region  and  matching  the 
tangential  components  of  E  and  H  at  the  interfaces  to  solve  for  the  unknown  constants. 
For  two  layers,  there  are  six  unknown  constants  to  be  determined  by  solving  a  sixth-order 
system  of  equations.  While  this  is  a  valid  means  of  proceeding,  the  derivation  is  lengthy 
and  difficult,  and  it  is  impractical  for  three  layers  or  more. 

For  a  plane  wave  polarized  either  TE  or  TM  with  respect  to  the  surface  normal  and 
propagating  along  the  line  defined  by  the  spherical  angles  {di,(pi),  the  field  anywhere  in  an 
,V-layer  dielectric  medium  can  be  derived  by  the  following  steps.  First,  rotate  the 
coordinate  axes  so  that  the  incident  propagation  vector  is  in  the  x-z  plane  of  the  new 
coordinate  system.  The  incident  angle  (pi  is  equal  to  the  (p  of  Fig.  3.2.  If  the  incident  field 
is  derived  from  the  unit  vector  potentials  A  and  F,  then  the  plane-wave  fields  propagating 
in  a  free-space  region  in  the  absence  of  scatterers  are  given  by  Harrington  [24], 

eTo^*  =  [  -x'jky  +  'yjk,-  ,  (3.54a) 

Er  =  7^1  i'jy  +  yjk,.  ]/*.'■ "  *  V  ,  (3,54b) 

Applying  the  coordinate  transformation  (Eq.  (3.2)), 

=  y  [  -jky-  sin^  cos(p  ~  ,  (3.55a) 

cos<;^  +  jky  sin (p  ~  ,  (3.55b) 


where 


kx  -  k^‘<^os<p  +  ky-sinp  , 


ky  =  0. 

As  before,  the  primed  variables  x\  y',  z’  refer  to  the  original  coordinate  system,  and  the 
unprimed  variables  x,  y,  z  refer  to  the  rotated  coordinate  system.  Second,  translate  the 
fields  problem  into  an  equivalent  transmission-line  problem  by  assigning  a  characteristic 
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impedance  and  propagation  constant  yi  to  each  layer,  where  the  polarization  p  =  TE  or 
TM. 


Vi  -  Eriko  , 

yTE  _  ~^y  _ 

~  7i  ’ 

7™  -  £l  = 

Hy  jCDE^iEo  • 

Third,  alternating  between  IJz)  and  Z(z)  from  the  load  end  of  the  transmission  line  to  the 
source  end,  calculate  the  generalized  reflection  coefficient  at  each  interface.  Remember  that 
although  V(z),  I{z),  and  Z(z)  are  continuous  at  the  interfaces,  r\z)  is  not.  For  the  two-layer 
example  introduced  in  Sec.  3.1, 

^02  +  ^TOCOth72^2 

~  '^02”T - - -  > 

^00  Zo2COth72^2 


zf„  =  Zgi 


Zqi  +  Z^cothyiti 

Z^  +  Z^icoth/iti 


z;+zSi 

yP  _  yP 

Zfn  +  ZSl 


yP  _  yP 


00 


zf„  +  Z^ 


00 


TE  TM 

Since  Ey  ^  -V  and  <->  +V,  the  complex  amplitude  of  the  incident  voltage  wave  Vq  is 
given  by 
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=  J^y-  sm<p  +  jk^>  cos</>  ,  (3.56a) 

^  ^  •  (3.56b) 

ycoeo 

Fourth,  compute  the  voltage  in  each  region  beginning  at  the  source  end  of  the  transmission 
line.  For  the  two-layer  example,  in  region  0, 


V^(-ri)  =  VgVo'i  [  1  +  rg(-ri)  ]  . 

In  region  1, 

v^K)f21 

‘  1  +  ' 

vP{Q)  =  vf  [  1  +  r^(o)  ] . 

Fifth,  make  the  association  from  voltage  back  to  the  electric  field, 

-v'^^(o) , 

£™  =  +V™(0), 

and  sixth,  rotate  back  to  the  original  coordinate  system; 


(3.57) 

(3.58) 

(3.59) 

(3.60a) 

(3.60b) 


E|"‘^(x',y')  =  [  x'sin^  -  y'cos^lE^o  \  TE  inc.,  (3.61a) 

=  [  x'cos0  +  y‘sin(^]E™/*^'*'^  TM  inc.  (3.61b) 


This  completes  the  derivation  of  the  incident  field. 
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CHAPTER  4 

THEORETICAL  AND  NUMERICAL  CONSIDERATIONS 

Before  the  equations  developed  in  the  previous  chapter  can  be  implemented  on  the 
computer,  there  are  a  number  of  points  that  must  be  considered. 


4.1  General  Characteristics  of  the  Inner  Products  Arising  in  the  Analysis 
of  Finite  Frequency  Selective  Surfaces 

If  a  function  iyoUj)  is  shifted  relative  to  the  origin  by  an  amount  (xj.yj),  where 
ijix.y)  }jQ(x-Xj,y-yj),  then  the  Fourier  transforms  of  the  two  functions  are  related  by  the 
translation  property  of  the  Fourier  transform.  Modifying  the  notation  slightly  so  that 
jjQ{a,P)  is  real,  the  Fourier  transform  of  jy(x,y)  can  be  expressed  in  the  form, 

ij(a,p)  =  KjijQia.p)  ,  (4.1) 

where  Kj  is  a  complex  constant.  If  the  original  function  was  either  real  and  even  or  real 
and  odd,  ]  jQ(a,p)  will  be  either  real  and  even  or  real  and  odd.  This  is  true  for  the  vast 
majority  of  basis  and  testing  functions  used  in  practice,  one  exception  being  the  traveling- 
wave  expansion  used  in  the  analysis  of  electrically-large  bodies. 

If  the  current  density  on  a  finite-width  two-dimensional  surface,  for  example,  is 
discretized  using  entire-domain  basis  functions,  such  that  yj  =  0  for  all  j,  then  substituting 
representation  (4.1)  into  Eq.  (3.47),  the  matrix  element  Zy  is  given  by 


-K^Kj 

iTtTy 


X  G{a,p„)  ijoioc.p^)  ei°^‘>da 


(4.2) 
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where  Xij  =  JCj  -  xj,  and  the  superscript  signifies  transpose.  Since  j  io((X,p)  is  real,  the 
complex  conjugate  appears  only  on  the  complex  constant  Ki  and  on  the  phase  shift  term 

If  (p(a)  is  defined  as  the  product  j.Q(a,p)G(a,Pn)3jQ((X,Pn)^  which  is  a  complex¬ 
valued,  scalar  function  of  a,  except  for  isolated  singularities  due  to  the  Green's  function 
(which  is  the  subject  of  the  next  section),  <p(a)  is  relatively  smooth.  The  term  on  the 
other  hand,  which  depends  on  the  separation  between  the  basis  and  testing  functions,  can 
be  highly-oscillatory,  and  in  general,  the  greater  the  distance  between  a  pair  of  basis  and 
testing  functions,  the  higher  the  oscillation.  A  plot  of  the  integrand  is  shown  in  Fig.  4.1 
for  a  finite  strip  grating  whose  current  density  flowing  perpendicularly  to  the  strips  has 
been  discretized  with  triangular  basis  functions  with  a  half  width  A  of  IVIO. 


(d,(p)  =  (.1,  .1  ) 


alpha  (1/m) 

Figure  4.1.  Plot  of  the  integrand  representing  the  mutual  coupling  between 
two  .x-directed  triangular  currents  separated  by  2Aq. 


If  the  rate  of  oscillation  is  high  enough  relative  to  the  variation  of  (pa),  in  the  integral,  each 
positive  half  cycle  of  the  integrand  will  cancel  an  adjacent  negative  half  cycle  of  the 
function,  making  the  result  zero,  or  nearly  zero.  This  intuitive  observation  is  stated 
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mathematically  in  the  Riemann-Lebesgue  Lemma  [13],  which  says,  if  (p{a)  is  absolutely 
integrable  on  an  interval  (a.b),  then 

b 

lim  f  (p(oc)  =  0  , 

Ja 

where  a  and  b  are  finite  or  infinite  constants.  It  is  not  necessary  for  (p{oc)  to  have  a 
bounded  derivative  for  the  above  result  to  hold. 

A  sufficient  condition  to  guarantee  that  (p{a)  is  absolutely  integrable  is  that  each 
dielectric  layer  be  slightly  lossy.  Otherwise,  the  integrand  will  contain  non-integrable 
singularities  corresponding  to  unattenuated  guided  modes  in  the  dielectric-slab  waveguide. 

The  following  conclusions  are  deducible  from  the  above  premises.  (1)  If  the 
dielectric  layers  have  some  loss,  however  small,  the  mutual  coupling  between  FSS 
elements  will  become  negligibly  small  as  their  separation  increases.  This  suggests,  (2)  that 
beyond  a  cenain  separation,  it  will  be  possible  to  neglect  the  direct  mutual  coupling 
between  two  elements  in  the  array.  In  very  large  arrays,  it  is  only  necessary  to  calculate  the 
mutual  coupling  between  one  element  and  a  finite  number  of  its  nearest  neighbors.  The 
remaining  terms  will  be  negligible  by  comparison  and  may  be  set  to  zero  without 
introducing  significant  error. 

The  coupling  distance,  or  maximum  separation  for  which  direct  mutual  coupling  must 
be  taken  into  account,  depends  on  the  parameters  of  the  dielectric  support.  If  the  dielectrics 
are  extremely  lossy,  for  example,  the  fields  produced  by  one  element  will  attenuate  rapidly 
away  from  the  source  and  will  have  little  direct  effect  on  the  currents  induced  on  a  distant 
element  of  the  array.  If  the  elements  are  part  of  a  free-standing  surface,  numerical 
experiments  have  shown  that  if  the  x  and  y  periodicities  are  roughly  equal  to  A/2, 
significant  mutual  coupling  exists  only  between  elements  that  are  less  than  about  seven 
periods  apart.  For  surfaces  supported  by  dielectric  layers,  or  backed  by  perfectly- 
conducting  ground  planes,  the  coupling  distance  will  be  somewhat  higher. 
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A  third  conclusion  deducible  from  the  above  premises  is  the  following:  for  distant 
interactions,  the  dominant  coupling  mechanism  is  the  guided  modes  in  the  dielectric  layers. 
U  Xij  is  large  but  not  infinite,  over  intervals  where  (p{a)  is  relatively  smooth,  as  already 
shown,  the  factor  will  cause  cancellation  to  occur  in  the  evaluation  of  the  integrals. 
For  large  distances  Xij,  therefore,  significant  contributions  to  the  inner  products  will  occur 
when  <p(a)  is  singular  or  sharply  peaked.  These  features  are  associated  with  the  guided 
modes  in  the  dielectric-slab  waveguide. 


4.2  Singularities  of  the  Spectral  Green's  Function 

The  analysis  of  periodic  structures  by  the  spectral-domain  method  calls  for  sampling 
the  product  of  the  Green's  function,  basis,  and  testing  functions,  in  the  Fourier-transform 
domain,  at  the  transverse  wave  numbers  oc^,  of  the  Floquet  space  harmonics.  These 
are  represented  graphically  by  the  intersections  of  the  grid  lines  in  the  reciprocal  lattice  of 
Fig.  4.2.  Since  the  frequency,  permittivities,  thicknesses,  periodicities,  and  incident  angles 
are  chosen  independently,  it  is  unlikely  that  the  Green’s  function  will  be  sampled  .n  one  of 
its  singular  points,  represented  by  the  concentric  circles  in  Fig.  4.2.  However,  in  the 
unlikely  event  that  this  does  happen,  any  of  the  above  parameters,  for  example  the  incident 
angles,  can  be  adjusted  by  a  small  amount  to  shift  either  the  locations  of  the  singularities  or 
the  locations  of  the  sample  points  to  avoid  sampling  directly  on  a  singularity  of 

the  Green's  function. 


For  a  frequency  selective  surface  of  finite  width,  the  formulation  calls  for  a  number  of 

“7  mt 

integrations  to  be  performed  along  the  diagonal  lines  ^  -  («  -  cotG  -t-  ky, 

labeled  /i  =  0,  ±1,  ±2,  ...  in  Fig.  4.2.  For  free-standing  surfaces,  the  free-space  Green's 
function  will  contribute  an  integrable  singularity  to  the  integrand  at  k^^  =  kg^.  In  this  case, 
the  function  can  be  integrated  with  a  judicious  choice  of  a  numerical  integration  routine.  If 


the  printed  surface  is  deposited  on  the  surface  of  one  or  more  lossless  dielectric  layers,  for 
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I  t  »  t 

m=-l  m  -0  m  =  \ 

Figure  4.2.  Reciprocal  lattice  for  periodic  surfaces. 


a  real,  at  least  one  integration  will  pass  through  a  non-integrable  singularity  of  the  Green's 
function.  Viewing  the  situation  in  the  complex-a  plane,  the  path  of  integration  is  along  the 
real-a  axis,  where  a  number  of  singularities  are  encountered  for  a  small  number  of 
integrations  around  n  =  0  (Fig.  4.3(a)). 

There  are  at  least  two  methods  of  handling  the  singularities  [25].  If  the  dielectric 
layers  are  assumed  to  be  slightly  lossy,  the  singularities  migrate  off  the  path  of  integration 
as  shown  in  Fig.  4.3(b),  and  there  is  no  longer  any  theoretical  difficulty  in  performing  the 
integrations.  This  is  not  usually  an  inconvenience,  since  all  physical  materials  have  some 
loss.  However,  the  lossless  situation  can  be  considered  in  the  limit  as  the  loss  tends  to 
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lm(  a  }  Im(  a  } 


(c)  (d) 

Figure  4.3.  Complex- a  plane: 


(a)  singularities  of  the  spectral  Green's  function  for  a  lossless 
dielectric,  assuming  one  guided  mode  is  supported  in  the 

dielectric  (time  convention: 

(b)  locations  of  the  singularities  assuming  the  dielectric  has  a  small 
amount  of  loss. 

(c)  contour  of  integration  obtained  by  considering  the  lossless  case 
as  the  zero-loss  limit  of  (b), 

(d)  rectilinear  path  obtained  by  deforming  the  contour  of  (c). 

zero,  as  long  as  the  contour  of  integration  is  defined  in  such  a  way  as  to  encircle  the 
singularities  in  a  manner  which  is  consistent  with  the  lossy  case  (Fig.  4.3  (c)).  However, 
since  the  locations  of  the  singularities  can  only  be  determined  numerically,  and  since  the 
contour  can  be  deformed  through  any  analytic  region  of  the  function,  it  is  more  convenient 
to  integrate  along  a  rectilinear  path  as  shown  in  Fig.  4.3  (d). 
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Both  of  these  methods  were  found  to  be  successful  means  of  handling  the 
singularities.  Contour  integration  may  offer  some  advantages  in  terms  of  efficiency.  If  the 
distance  5  from  the  path  of  integration  to  the  singularities  is  varied  as  a  parameter,  then  the 
greater  the  distance,  the  smoother  the  integrand,  and  the  fewer  number  of  points  are 
required  to  integrate  the  function  to  the  same  accuracy.  In  Table  4.1,  a  typical  integral 
ansing  in  the  study  of  the  finite  strip  grating  is  evaluated  by  the  contour  integration  method, 
for  various  values  of  5.  The  current  density  flowing  parallel  to  the  axis  of  the  strips  has 
been  discretized  using  pulse  basis  functions.  The  dominant  self-interaction  term  is 
evaluated  using  an  adaptive  numerical  integration  routine  that  automatically  adjusts  the 
sampling  rate  in  order  to  accommodate  for  changes  in  the  smoothness  of  the  function. 
Table  4.1  shows  the  number  of  function  evaluations  required  to  achieve  the  same  reladve 
accuracy  of  0. 1  per  cent.  As  indicated  in  the  table,  the  least  number  of  function  evaluations 
is  exactly  one  third  of  the  greatest  number. 


Table  4.1.  Number  of  function  evaluations  required  to  integrate  a  function 
to  the  same  relative  accuracy  using  the  contour  integration  method. 
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The  results  of  Table  4.1  suggest  that  it  may  be  possible  to  find  an  optimum  path 
which  would  require  a  minimum  of  computational  effort.  Kong  and  Shen  [26],  to  find  the 
fields  of  an  antenna  submerged  in  a  dissipative  dielectric  medium,  define  vertical  branch 
cuts  which  are  actually  the  steepest-descent  paths  for  a  particular  configuration  of 
transmitting  and  receiving  antennas.  Since  the  integrand  is  exponentially  decaying  in  the 
direction  perpendicular  to  the  real- ct  axis,  it  can  be  integrated  using  only  a  few  evaluations 
of  the  function. 

The  finite  FSS  will  require  considerably  more  integration  time,  but  the  considerations 
are  identical  to  those  in  this  case. 

4.3  Numerical  Integration 

As  was  shown  in  previous  chapters,  the  spectral-domain  formulation  of  the  method 
of  moments  applied  to  the  analysis  of  finite  frequency- selective  surfaces  leads  to  a  number 
of  integrals  involving  the  Fourier  transforms  of  the  Green's  function,  basis,  and  testing 
functions.  Because  of  the  complexity  of  the  integrands,  the  integrals  cannot  be  evaluated 
exactly  and  some  method  of  numerical  integration  must  be  employed. 

An  intelligent  choice  of  a  numerical  integration  scheme  depends  on  the  nature  of  the 
integrand.  Integrals  with  special  difficulties  such  as  infinite  ranges  of  integration, 
singularities,  or  highly-oscillatory  integrands  should  be  treated  in  special  ways,  developing 
specific  automatic  routines  tailored  for  each  type  of  difficulty.  “The  programmer,” 
according  to  Davis  and  Rabinowitz,  is  “in  the  position  of  the  sportsman  who  has  given  up 
worms  but  now  must  know  which  fly  to  use  to  catch  his  frsh”  [27].  A  poor  choice  could 
lead  to  excessive  computation  time  or  give  erroneous  results. 

The  integrands  encountered  in  the  spectral-domain  method  of  moments  have  all  of  the 
above  difficulties  associated  with  them.  First,  in  the  lossless  limit,  for  <  Emiax 
integrands  have  a  number  of  isolated  singularities.  There  are  branch  points  where  the 
phases  of  the  propagation  constants  Yi  arc  discontinuous  and  singularities  caused  by  zeros 
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in  the  denominators  of  the  input  impedances  and  Even  if  the  singularities  are 

avoided  either  by  introducing  a  small  loss  into  each  region  or  by  deforming  the  path  of 
integration  into  the  complex  plane,  the  integrands  will  still  vary  rapidly  in  the  vicinities  of 
the  singularities.  Second,  interactions  between  basis  and  testing  functions  which  are 
separated  by  a  few  periods  are  characterized  by  highly-oscillatory  integrands.  As  the 
distance  between  a  pair  of  basis  and  te:  ting  functions  increases,  the  rate  of  oscillation 
increases  proportionately.  Third,  the  limits  of  integration  are  infinite. 

Of  the  three,  the  singularities  are  considered  to  be  the  most  serious  difficulty.  The 
oscillation  is  manageable  for  arrays  with  moderate  dimensions,  while  for  large  arrays,  the 
integrations  can  be  confined  to  the  regions  around  the  poles.  If  the  integrand  has  a  certain 
minimum  asymptotic  rate  of  decay,  the  limits  of  integration  can  be  truncated  to  finite  values 
and  the  integral  evaluated  over  the  finite  interval. 

The  sampling  rate  in  numerical  integration  is  determined  by  the  number  of  points 
necessary  to  represent  the  integrand  in  a  given  region.  For  kp  <  £rmax  .  ^he  sampling 

rate  is  dominated  by  the  presence  of  the  singularities.  Outside  this  region,  the  oscillatory 
behavior  of  the  integrand  dominates.  Because  the  integrand  has  two  regions  dominated  by 
different  behaviors,  because  the  integrand  contains  a  number  of  singularities,  and  because 
the  locations  of  these  singularities  can  only  be  determined  numerically,  an  adaptive 
numerical  integration  routine  that  automatically  adjusts  the  sampling  rate  to  accommodate 
for  changes  in  the  smoothness  of  the  integrand  is  a  natural  choice. 

The  advantage  of  such  a  scheme  is  that  the  integrand  is  sampled  heavily  only  in 
regions  where  the  function  is  varying  most  rapidly.  If,  over  a  given  interval,  the  integrand 
is  relatively  smooth,  the  sampling  rate  is  automatically  reduced  to  a  minimum.  The 
algorithm  is  also  capable  of  handling  singularities.  Davis  and  Rabinowitz  report  that  an 
adaptive  routine  based  on  Simpson's  rule  has  been  used  to  successfully  integrate  1/vTri 
through  zero.  Hillstrom  [28]  has  done  some  numerical  experiments  that  indicate  that  if  the 
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integrand  is  sharply  peaked,  or  if  great  accuracy  is  desired,  a  high-order  adaptive  routine  is 
most  efficient. 

Many  variations  of  adaptive  integration  have  been  proposed  in  the  literature.  Rice, 
who  is  quoted  by  Davis  and  Rabinowitz,  has  studied  the  subject  in  detail  and  estimates  that 
there  are  about  a  million  adaptive  quadrature  algorithms  which  are  potentially  reasonable 
and  interesting.  While  a  complete  discussion  of  the  subject  would  be  out  of  place  here,  the 
essential  features  of  adaptive  integration  are  embodied  in  the  following  example. 

On  a  given  interval,  the  function  is  first  integrated  using  a  previously-chosen 
fundamental  integration  rule.  The  interval  is  then  halved  and  the  same  rule  is  applied  to 
each  of  the  two  subintervals.  If  the  result  of  integrating  over  the  two  subintervals  is  within 
a  specified  tolerance  of  the  result  for  the  total  interval,  the  result  of  the  two  integrations  is 
accepted,  and  a  new  interval  is  selected  for  the  next  iteration.  If  the  result  of  the  two 
integrations  does  not  fall  within  the  allowable  tolerance,  one  of  the  two  subintervals  is 
selected  as  the  total  interval  for  the  next  iteration,  and  the  process  is  repeated.  For  example, 
if  ,VCn  is  the  result  of  applying  an  n-point  Newton-Cotes  integration  rule  over  [a, 6],  and 
2  X  NCn  is  the  result  of  applying  the  same  rule  to  each  of  the  two  subintervals,  then,  if 

I  2  X  NC^-NC„  I  <  e I  2  X NC„  1 ,  (4.3) 

the  result  2  x  NCn  is  accepted  as  the  value  of  the  integral  over  [a,b].  Otherwise,  each  of  the 
subintervals  is  treated  as  the  original  interval  using  £/V2  as  the  pro-rated  error  tolerance, 
since  e/2  is  usually  too  strict  in  practice.  The  actual  error  is  usually  much  less  than  e. 

The  choice  of  an  integration  rule  to  be  used  as  the  basis  for  an  adaptive  algorithm 
depends  on  a  number  of  considerations.  The  theory  of  approximate  quadratures  depends 
on  replacing  the  integral  \f{x)dx  by  another  integral  (p(x)dx,  where  q){x)  can  be 
determined  in  a  simple  way  [29].  The  function  (p{x)  is  expressed  as  a  series  of  weighted 
orthogonal  functions  over  [a,b]',  then  (p(x)  is  integrated  exactly  to  obtain  an  approximation 
to  the  desired  integral.  If  the  integrand  is  sampled  at  the  (n  +  1)  points  yo>  Ji.  Vn,  the 
result  of  the  integration  can  be  expressed  in  the  form 
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/  (pix)  dx  =  /4oyo  +  ^  iJi  +  .  .  .  +  ,  (4.4) 

•'a 

where  the  («  +  1)  constants  Af^  are  independent  of  the  (n  +  1)  abscissas  y^-  In  the 
Newton-Cotes-type  formulas,  the  points  ym  are  fixed.  It  follows  that  if  fix)  is  a 
polynomial  of  degree  <  n,  the  error  in  integrating  (p(x)  can  be  made  to  vanish  by  the  proper 
choice  of  the  /4^'s.  In  Gauss's  method,  both  the  sample  points  y^  and  weights  are 
determined  in  such  a  way  that  the  difference  between  ^f{x)dx  and  |  (p{x)dx  is  a  minimum. 
Since  there  are  now  (In  +  2)  available  constants,  an  (n  +  l)-point  Gaussian-quadrature 
formula  is  capable  of  integrating  exactly  a  polynomial  of  degree  <2n  +  1.  Thus,  for  the 
same  relative  accuracy,  a  Gaussian-quadrature  formula  requires  fewer  points,  or,  for  a 
given  number  of  sample  points,  a  Gaussian-quadrature  formula  is  more  nearly  exact  than 
the  corresponding  Newton-Cotes  formula. 

Integration  rules  based  on  equally  spaced  abscissas,  such  as  the  Newton-Cotes-type 
formulas,  have  the  advantage  that  as  each  interval  is  halved,  only  half  of  the  function  values 
need  to  be  computed,  since  the  other  half  is  available  from  a  previous  iteration.  Many  high- 
order  Newton-Cotes  rules  contain  weights  of  opposite  signs,  which,  because  of  finite  word 
length,  could  result  in  a  loss  of  accuracy.  Newton-Cotes  rules  based  on  an  even  number  of 
points  have  rarely  been  used  in  practice,  since  they  are  generally  less  accurate  than  the  next 
lower  order  rule  of  the  same  type,  which  can  be  verified  by  examining  the  error 
expressions  given  by  Abramowitz  and  Stegun  [30].  A  notable  exception  to  this  is  the 
trapezoidal  mle. 

As  stated  before,  Hillstrom  has  shown,  by  numerical  experiment,  that  if  the  integrand 
is  sharply  peaked,  or  if  great  accuracy  is  desired,  an  adaptive  high-order  rule  is  most 
efficient.  Hillstrom's  experiments  were  repeated  using  the  adaptive  procedure  outlined  in 
this  section.  Adaptive  algorithms  ANC3,  ANC7,  and  ANCl  1,  based  on  3,  7,  and  1 1-point 
Newton-Cotes  formulas,  were  applied  to  two  different  test  integrals  having  sharply  peaked 
or  singular  integrands.  ANC3  is  an  adaptive  Simpson's  rule,  which  is  the  minimum- 
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degree  rule  used  adaptively  in  practice.  ANC7  is  based  on  a  7-point  Newton-Cotes  rule, 
which  is  the  highest-degree  rule  based  on  an  odd  number  of  points  with  all  positive 
weights,  and  ANCl  1  is  based  on  the  highest-degree  rule  given  by  Abramowitz  and  Stegun 
[30], 

Each  of  the  routines  was  required  to  integrate  the  same  function  to  the  same  relative 
accuracy  £,  and  each  was  limited  to  the  same  minimum  sampling  interval,  which  was 
chosen  to  be  approximately  2/10®.  If  convergence  had  not  been  achieved  at  this  point,  the 
most  accurate  estimate  was  accepted  as  the  value  of  the  integral  over  that  subinterval  and  the 
process  was  allowed  to  continue.  An  algorithm  is  considered  most  efficient  if  it  requires 
the  least  number  of  functional  evaluations  to  obtain  the  result. 

The  first  integral  is  of  the  form  J.\(ac2  +  p2)-i  where  is  a  parameter.  For  small 
p.  the  integrand  has  a  peak  of  height  at  the  origin  and  is  approximately  equal  to  one  at 
the  endpoints.  Table  4.2  shows  the  number  of  function  evaluations  required  by  each 
routine  to  achieve  the  same  relative  accuracy.  Results  are  also  provided  for  the  popular 
ROMBERG  integration  method,  to  highlight  the  advantages  of  using  adaptive  integration 
for  sharply  peaked  integrands. 

Although  the  adaptive  algorithms  described  here,  and  in  particular  the  convergence 
criterion,  are  somewhat  different  from  those  implemented  by  Hillstrom,  the  results  are  very 
similar.  The  “best”  integration  rule  to  use  in  an  adaptive  scheme  depends  both  on  the 
desired  accuracy  and  on  the  smoothness  of  the  integrand.  If  extreme  accuracy  is  desired, 
an  adaptive  routine  built  on  a  high-order  rule  is  most  efficient;  but  if  two  or  three  digits  of 
accuracy  is  sufficient,  for  the  range  of  the  parameters  tested,  an  adaptive  routine  built  on  a 
moderately  high-order  Newton-Cotes  formula  ANC7  is  most  efficient.  The  comparison 
with  the  results  of  the  ROMBERG  algorithm  is  striking. 


Table  4.2.  Comparison  of  several  adaptive  quadrature  routines  for 
integrating  sharply-peaked  functions  of  the  form  {x^  -t- 
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A  second  type  of  difficulty  encountered  in  the  analysis  of  finite  frequency  selective 
surfaces  is  an  integrable  singularity  of  the  type  l/Vx^  -  1.  This  type  of  singularity  is 
present  in  the  transform-domain  expression  for  the  free-space  Green's  function,  Eqs.  (3.3) 
and  (3.26).  At  the  singular  points,  x  =  ±1,  the  magnitude  of  the  integrand  becomes 
unbounded,  while  the  function  reverses  from  real  to  imaginary  and  vice  versa.  This 
function  can  be  integrated  exactly  to  give 


This  integral  was  evaluated  using  the  same  adaptive  quadrature  schemes  as  before.  The 
results  are  summarized  in  Table  4.3.  Of  the  routines  tested,  for  0.01  <  e  <  0.005,  the 
adaptive  Simpson's  rule  was  most  efficient;  for  0.002  <  e  <  0.0001,  the  adaptive  7-point 
rule  was  most  efficient.  As  in  the  first  example,  for  extremely  tight  error  bounds,  say£  < 

10- ^,  the  higher-order  rule  ANCll  is  expected  to  become  more  efficient  than  the  lower- 
order  mles  ANC3  and  ANC7. 

These  conclusions  are  based  on  the  performance  of  integration  rules  with  evenly- 
spaced  abscissas.  The  question  remains  as  to  whether  there  is  any  advantage  in  using  an 
integration  rule  with  non-equally-spaced  abscissas,  such  as  the  Gauss-type  formulas. 
Chugh  and  Shafai,  for  example,  report  that  an  iterative  nonadaptive  algorithm  built  on  an 

1 1 - point  Gauss-Legendre  integration  rule  was  more  accurate  and  more  efficient  than  the 
popular  ROMBERG  integration  method  at  handling  a  particular  two-dimensional  phase 
integral  arising  in  electromagnetics  [31]. 

The  fact  that  the  abscissas  are  unequally-spaced  makes  an  adaptive  Gaussian- 
quadrature  algorithm  particularly  robust.  When  a  given  interval  is  bisected,  since  the 
abscissas  are,  for  example,  given  by  the  zeros  of  Legendre  polynomials,  none  of  the 
sample  points  coincide  with  previous  values,  and  subsequent  integrals  are  based  on  entirely 
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Table  4.3.  Comparison  of  several  adaptive  quadrature  routines  for 
integrating  a  singularity  of  the  form  -  1  . 
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different  values  of  the  function.  If  the  convergence  criterion  is  satisfied  under  these 
conditions,  the  result  is  especially  reliable. 

However,  the  source  of  all  the  advantages  of  using  a  Gaussian-quadrature  formula  in 
an  adaptive  routine  is  also  the  cause  of  its  main  disadvantage.  Since,  in  general,  no  two 
abscissas  coincide,  all  previous  function  evaluations  are  discarded  in  the  next  stage  of 
refinement.  The  overall  efficiency  of  an  adaptive  Gaussian-quadrature  routine,  relative  to  a 
comparable  Newton-Cotes  routine,  will  therefore  depend  on  the  relative  importance  of  the 
following  considerations:  (1)  the  number  of  points  used  in  actually  evaluating  the  integral 
is  expected  to  be  somewhat  less  for  an  adaptive  Gaussian-quadrature  routine,  but  (2)  the 


cost  of  determining  these  points,  which  is  equal  to  the  cost  of  discarding  all  previous 
function  values,  may  be  prohibitive.  Clearly,  the  second  factor  depends  on  the  complexity 
of  the  particular  integrand  being  considered.  As  the  integrand  becomes  more  difficult  to 
compute,  the  added  computation  time  due  to  calculating  and  then  discarding  previous 
function  values  will  eventually  outweigh,  then  overwhelm,  the  time  saved  in  using  a  more 
accurate  integration  rule.  Unfortunately,  there  is  no  way  to  determine  the  complexity  of  the 
integrand  where  this  transition  occurs. 

A  numerical  experiment  was  conducted  in  order  to  determine  the  efficiency  of  an 
adaptive  Gaussian-quadrature  routine  relative  to  a  comparable  Newton-Cotes  routine  for  a 
particular  integral  arising  in  the  analysis  of  the  finite  strip  grating.  Since  Gaussian- 
quadrature  formulas  are  open  type,  i.e.,  since  the  function  is  not  evaluated  at  the  endpoints 
of  the  interval,  a  10-point  Gauss-Legendre  rule  was  compared  to  an  11 -point  Newton- 
Cotes  rule.  In  order  to  integrate  a  function  of  this  complexity,  the  adaptive  Gaussian- 
quadrature  routine  required  from  two  to  three  times  as  long  as  the  comparable  Newton- 
Cotes  algorithm. 

4.4  Symmetry 

The  input  impedances  and  that  are  used  in  the  derivation  of  the  spectral- 
domain  Green's  function  are  functions  of  +  ^92  only,  and  are,  therefore,  even  functions 
of  both  variables.  According  to  Eqs.  (3.19)-(3.21),  Gxx  and  Gyy  are,  therefore,  even 
functions  of  ct  anu  n,  vvliiie  Gx^  is  ?<i  odd  luacticn  of  both  variables.  Since  nearly  all  of 
the  basis  and  testing  functions  used  in  practice  are  either  even  or  odd,  their  Fourier 
transforms  are  either  even  or  odd  as  well.  In  this  section,  we  will  investigate  under  what 
conditions  symmetry  can  be  used  to  reduce  a  two-dimensional  summation  or  integration 
over  the  infinite  a-P  plane  to  a  summation  or  integration  over  the  half  plane  a  >  0  or  the 
quarter  plane  a  >  0,  >  0,  reducing  the  total  number  computations  to  one  half  or  one 

quarter  of  its  original  value. 
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To  take  advantage  of  symmetry,  the  integrand  or  summand  must  be  evaluated  at 
symmetrical  points  (±a,  ±/3).  Referring  to  Fig.  4.2,  it  is  clear  that  for  the  periodic  FSS, 
the  sample  points  of  the  summand  are  not  symmetrical,  in  general,  about  either  axis,  and 
symmetry  cannot  be  used  except  in  special  cases.  If  the  skew  angle  Q  is  equal  to  90°  and 
the  incident  angle  (p  is  zero,  then  the  sample  points  of  G,  j  i,  and  j  n  =  0,  ±1, 

±2, ...,  will  be  symmetrical  with  respect  to  the  a  axis,  and  the  computation  time  in  filling 
the  moment  matrix  would  be  reduced  by  approximately  one  half.  If,  in  addition,  0  =  0, 

^TDTl 

then  the  sample  points  ,  /n  =  0,  ±1,  ±2, ...,  would  also  be  symmetrical  about  the 

P  axis,  allowing  the  number  of  computations  to  be  further  reduced  to  approximately  one 
fourth  the  original  number.  However,  if  the  incident  field  is  derived  from  the  magnetic  and 
electric  vector  potentials,  Eqs.  (2,l)-(2.2),  the  field  for  (d,(p)  =  (0,0)  is  not  defined. 

On  the  other  hand,  if  the  problem  of  interest  is  a  frequency  selective  surface  of  finite 
width,  rather  than  summations  at  the  discrete  points  {cc^,pfnn),  the  formulation  calls  for 
integrations  along  the  diagonal  lines  /?„  =  -  («  *  kx)cotQ  +  ky,  shown  in  Fig.  4.2. 

As  in  the  periodic  case,  generally  the  integrands  are  not  evaluated  at  symmetrical  points  and 
symmetry  cannot  be  used.  However,  there  are  numerous  practical  situations,  less 


restrictive  than  in  the  periodic  case,  where  symmetry  can  be  used.  If  12  =  90°,  then  the 

'2,7rfi 

lines  of  integration,  +  ICy,  are  horizontal,  and  the  even/odd  character  of  the 

integrands  with  respect  to  a  can  be  used  to  decrease  computation  time  by  one  half.  This  is 


true  for  the  strip  grating,  for  example,  and  for  a  large  number  of  other  geometries  used  in 


practice.  If,  in  addition,  (p  =  0,  then  ky  =  0,  and  the  lines  of  integration  will  be  located 
symmetrically  about  the  a  axis,  allowing  another  factor  of  two  reduction  in  computation 


time. 


For  the  finite-FSS  problem,  which  requires  a  full  two-dimensional  integration  over 
the  infinte  a-/3  plane,  if  the  basis  and  testing  functions  are  real  and  either  even  or  odd, 
symmetry  can  always  be  used  to  reduce  the  region  of  integration  to  the  quaner  plane  a  >  0, 
/3>0. 
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4.5  Translation,  Reciprocity,  and  Parallel  Matrix  Fill 

The  analysis  of  finite  frequency  selective  surfaces  by  the  spectral-domain  method  of 
moments  requires  an  enormous  amount  of  computation  time  unless  careful  consideration  is 
given  to  the  efficient  computation  of  the  matrix  elements.  It  has  already  been  shown  that  a 
judicious  choice  of  a  numerical  integration  technique  can  mean  an  order  of  magnitude 
speedup  in  generating  the  moment  matrix.  Equally  significant  is  the  computation  time 
saved  by  eliminating  redundant  calculations.  If  the  current  density  on  a  finite- width  FSS  is 
discretized  with  entire-domain  basis  functions  whose  unshifted  versions  are  real  and  either 
even  or  odd,  then  the  matrix  element  Zy  can  be  expressed  in  the  form 

oe 

^  n  =  -oo 

where  AT,-  and  Kj  are  complex  constants,  j,o  and  jyo  are  real,  and  (x/,  yi)  and  {Xj,  yp  are  the 
coordinates  of  the  basis  and  testing  functions  relative  to  the  origin.  According  to  this 
equation,  the  interaction  between  basis  function  j at  position  (xy,  yp  and  testing  function 
j  ,0  at  position  (x,-,  yp  depends  only  on  the  distances  x,y  =  (x,-  -  Xy  )  and  y,y  =  (y,-  -  yp.  For 
example,  if  the  current  density  on  a  finite-width  frequency  selective  surface  of  thin  linear 
dipoles  is  represented  in  terms  of  rooftop  basis  functions  as  shown  in  Fig.  4.4,  and 
Galerkin's  testing  procedure  is  used,  many  of  the  interactions  are  identical:  =  Z22,  etc., 

^12  ~  ^23’  ^13  ~  ^24’  ^16  ~  ^27>  Z17  =  Z2g,  CtC.,  Z|g  =  Z29,  CtC. 

If  the  basis  functions  are  numbered  as  in  Fig.  4.4  and  arranged  in  a  vector,  the 
resulting  matrix  Z  will  have  two  levels  of  Toeplitz  structure,  as  indicated  in  Fig.  4.5.  The 
full  matrix  will  be  a  block-Toeplitz  matrix,  and  each  Toeplitz  block  will  itself  be  a  Toeplitz 
matrix,  i.e.,  the  independent  elements  of  Z  lie  in  the  first  row  and  column  of  the  blocks 
and  the  remaining  elements  are  obtained  by  copying  these  elements  parallel  to  the  main 
diagonal.  The  result  is  that  instead  of  computing  (A/g)^  matrix  elements,  where  N  is  the 
number  of  elements  in  a  unit  row  and  /g  is  the  number  of  basis  functions  per  element,  it  is 
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Figure  4.4.  Rooftop  approximation  for  the  current  on  a  finite-width  array 
of  thin  dipoles. 


Figure  4.5.  Independent  elements  of  the  moment  matrix  Z  for  the  example 
of  Fig.  4.4,  showing  two  levels  of  Toeplitz  structure. 


only  necessary  to  compute  the  {2N  -  1)(2Iq  -  1)  independent  interactions  in  the  first  row 
and  column  of  each  Toeplitz  block.  The  actual  amount  of  computation  time  saved  will 
depend  on  the  size  of  the.  matrix. 
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The  translation  property  that  produced  the  significant  savings  in  this  example  is  not 
dependent  on  the  details  of  the  example.  The  following  general  conclusion  is  valid:  if  a 
finite  frequency  selective  surface  consists  of  a  uniform  array  of  identical  elements,  and  the 
unknown  on  each  element  is  discretized  with  the  same  set  of  basis  functions,  which  are 
ordered  in  the  same  way  and  grouped  according  to  element,  the  matrix  Z  will  be  a  block- 
Toeplitz  matnx.  The  second  level  of  Toeplitz  structure  in  this  example  is  due  to  the  choice 
of  subdomain  basis  functions;  however,  this  type  of  structure  can  also  arise  with  entire- 
domain  basis  functions  in  the  analysis  of  frequency  selective  surfaces  which  are  finite  in 
two  dimensions. 

A  different  matrix  structure  and  a  different  type  of  redundancy  show  up  in  the 
moment  matrix  when  the  current  on  a  finite  frequency  selective  surface  is  represented  in 
terms  of  entire-domain  basis  functions  with  Galerkin  testing.  It  is  easily  verified  that  the 
product  J/qG  jyo  is  independent  of  the  interchange  of  j  ,o  and  }  jq,  i.e.. 


This  suggests  that  the  inner  product  is  unchanged  if  the  functional  forms  of  the  basis  and 
testing  functions  are  interchanged,  provided  the  source  patch  and  test  patch  remain  the 
same.  For  example,  if  the  current  on  a  finite  array  of  thin  horizontal  dipoles  is  represented 
in  terms  of  a  trigonometric  series  as  shown  in  Fig.  4.6,  then  Z^2  -  ^21  >  ^13  =  ^31  >  ^23  = 
Z32,  etc.,  and  Z15  =  Z24,  Z^^  =  Z34,  Z26  =  Z35,  etc.;  i.e.,  each  Toeplitz  block  is 
symmetrical  with  respect  to  its  diagonal.  This  property  is  illustrated  in  Fig.  4.7. 


Figure  4.6.  Trigonometric  basis  functions  used  in  the  expansion  of  the 
current  on  a  three-element  frequency  selective  surface  of  finite 
width. 
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-ZHZ12Z13- 

Z14Z15Z16 

^17^18^19 

•  ^22^23 

■  ^25^26 

•  Z28Z29 

-  '  ■  ^32- 

-  ■  ■  ■^36- 

-  ■  ■  ■^39- 

Z41Z42Z43 
•  252Z53 

. 

-  ■  ^63- 

_  ■  •  •  _ 

*  *  '  — 

Z71Z72Z73 
■  ^82^83 

. 

-  '  ■  ^93- 

_  • 

_  ■ 

Figure  4.7.  Matrix  structure  for  the  example  of  Fig.  4.6,  showing  only  the 
unique  elements  of  the  matrix  (see  text). 


Under  certain  conditions,  it  is  possible  to  make  use  of  reciprocity.  If  the  skew  angle  = 
90°,  such  that  the  integrations  in  the  previous  example  are  performed  along  the  horizontal 
lines  /3n  =  —  +  ky,  then  for  the  even-even  interactions,  Zj3,  Z14,  Z15,  etc.,  AT,-  and  Kj  are 
real,  and  since  Gxx  is  an  even  function  of  a,  then 


Z..  = 


^  X  G„{a,l5„)  hjoia.Pr)  (4.6a) 


y  2nT. 


y  n  =  ~ 


iTtr,  ^ 


n  —  —00 


Jxio  (^xx  Jxjo  cosoxy  da 


(4.6b) 


and,  therefore,  Zjj  =  Zjj.  In  this  particular  example,  this  is  true  regardless  of  the  particular 
combination  of  basis  and  testing  functions  considered.  Therefore,  in  addition  to  the 
previous  reductions,  it  is  only  necessary  to  compute  the  elements  above  or  below  the  main 
diagonal.  The  resultant  matrix  is  shown  in  Fig.  4.8. 

In  addition  to  recognizing  that  many  of  the  matrix  elements  are  either  identical  or 
otherwise  simply  related,  large  reductions  in  computation  time  can  be  achieved  by 
computing  several  matrix  elements  in  parallel.  Referring  again  to  the  expression  for  the 
matrix  elements  for  a  finite-width  frequency  selective  surface. 
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^1 1^12^13 

^142^15^16 

^17^182^19 

■  ■^22^23 

■  ^25^26 

•  ^2S^29 

_  •  •  ■^33- 

.  ■  ■  Z36_ 

-  ■  ■  ^39- 

_  • 

Figure  4.8.  Matrix  structure  for  the  example  of  Fig.  4.6,  showing  the 
unique  elements  of  the  matrix  for  Q  =  90°  (see  text). 


^  2^'^'  21  G{a.(5„)  }jo{a,p^)  da  , 

^  n  =  -<»» 

the  idea  is  to  consider  the  product  j,oGjyo  ^  core  integrand  and  compute  all  of  the 
matrix  elements  at  once  for  various  values  of  (xi  -  xj)  and  (y,-  -  yj).  The  evaluation  of  G  at 
a  single  point  ia,p)  requires  solving  two  A'  layer  transmission-line  problems  (one  TE  and 
one  TM),  and,  depending  on  the  structure  of  the  program,  calculation  of  j  /  and  jy  could 
require  several  nested  calls  to  several  different  function  subprograms.  By  comparison,  the 
calculation  of  requires  very  little  time.  This  strategy  was  described  by 

Pozar  [32]  in  a  recent  paper  on  spectral-domain  analysis  of  a  fmite  microstrip-patch  array  of 
thin  conducting  dipoles. 

Unlike  when  the  method  of  moments  is  applied  in  the  space  domain,  the  interactions 
between  the  most  widely  separated  pairs  of  basis  and  testing  functions  are  the  most  difficult 
to  compute.  This  is  due  to  the  presence  of  the  highly-oscillatory  factor  ^j)'^JP(yi  ■  yj) 
in  the  integrand,  whose  rate  of  oscillation  increases  as  the  distances  {xi  -  xj)  and  iyi  -  yj) 
increase.  For  the  dominant  self-terms  of  the  matrix,  this  oscillatory  phase  factor  is  equal  to 
one;  therefore,  calculating  these  terms  requires  the  least  amount  of  computational  effort.  If 
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the  integrations  are  done  in  parallel,  such  that  all  of  the  integrands  are  sampled  more  heavily 
until  each  result  satisfies  the  exit  criterion,  then  the  dominant  self-terms  will  be  the  most 
accurate. 

Another  possible  meiuis  of  reducing  the  computation  time  to  fill  the  matrix  is  based  on 
the  observation  that  the  integral  is  of  the  form  of  an  inverse  Fourier  transform.  If  the 
separations  .t/,  are  all  integral  multiples  of  some  value,  then  the  FFT  can  be  used  to 
compute  all  the  interactions  for  a  given  pair  of  basis  and  testing  functions  simultaneously. 
Since  the  FFT  requires  that  the  samples  of  the  integrand  be  evenly  spaced,  this  is  equivalent 
to  ignonng  the  singularity.  This  strategy  has  been  used  by  Ko  and  Mittra  [33],  For  free¬ 
standing  surfaces,  the  error  involved  in  ignoring  the  singularity  does  not  appear  to  be  too 
great. 
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CHAPTER  5 
NUMERICAL  RESULTS 


5.1  Reflection  Coefficient  for  Periodic  Surfaces 

The  parameter  most  widely  used  to  characterize  the  scattering  response  of  frequency 
selective  surfaces  is  the  reflection  coefficient,  which  is  a  measure  of  the  reflectivity  of  the 
surface.  There  are  a  few  different  ways  of  defining  the  reflection  coefficient;  therefore,  it 
is  necessary  to  be  very  precise  in  the  definition  to  avoid  confusion. 

If  the  generalized  periodic  surface  of  Fig.  5.1  is  illuminated  by  a  plane  wave,  whose 
magnetic  or  electric  vector  potential  is  given  by 


then  the  scattered  electric  field  will  be  a  discrete  spectrum  of  plane  waves  of  the  form, 


E^(z  >z,)=  ^  +  (5-1) 

p.q  =  -«~ 

where  is  the  Kroniker  delta  function.  In  this  form,  the  scattered  field  can  be 

recognized  as  being  composed  of  two  parte:  (1)  the  scattered  field  associated  with  the 

vector  constants  is  due  to  the  induced  current  on  the  printed  surface  radiating  in  the 

presence  of  the  dielect  ic  medium;  and  (2)  the  scattered  field  associated  with  is  the  field 

uu 

reflected  from  the  dielectric  medium  with  the  printed  surface  removed. 

The  field  at  the  air-dielectric  interface  z  =  z,  and  the  induced  current  are  related  by  a 
spectral-domain  Green's  function,  similar  to  the  one  derived  in  formulating  the  integral 


equation, 


6? 


E^iz  =  z,) 

^xy 

4(z  =  0) 

Ey{z  =  zp 

1 

- \ 

>> 

1 - 

N 

II 

o 

1 _ 

(5.2) 


where  J  is  now  known  (or  an  approximation  to  J)  as  a  result  of  solving  the  discretized 

form  of  the  integral  equation.  The  derivation  of  the  spectral  Green's  function  is  similar  to 

the  derivation  in  Sec.  3.1,  and  the  details  are  omitted  here. 

z 


Figure  5.1.  Arbitrary  periodic  surface  surrounded  by  an  A-layer  dielectric 
medium. 


To  find  the  constants  the  strategy  is  to  take  the  Fourier  transform  of  Eq.  (5.1) 
and  match  the  two  representations  (5.1)  and  (5.2)  at  z  =  z,.  Recall  that  the  current  density  J 
is  related  to  tlie  current  density  on  a  single  element  Jq  by  the  following  relationship 


J(a,)3)  =  Jo(a-)3) 


Ik  In  'Y 
r„  sini2  ^ 


d{a  -  Upq)  5{p  -  pp^). 


(5.3) 


where  Upq  and  Ppq  are  given  by  Eqs.  (2.14).  Performing  the  indicated  operations  leads  to 

oo 

s  «  ~  -  Pp^) 

p.q  =  -^ 

oo 

Z-  -  9ir  9;r 

G^Clpq>Ppq)Ucip^.Ppq  >  T~'T~sinh  "  «pq^^^P  "  ^pq^^ 

p,q  =  -oo  ni  n  2 


or,  equating  the  two  series  term  for  term, 
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^pqt  ~  ^  >  (5.4) 

where  Cp^t  is  a  vector  of  the  transverse  (i.e.,  x-y)  components  of 

To  separate  the  scattered  field  into  its  TE  and  TM  components,  the  field  is  represented 
in  terms  of  z-directed  electric  and  magnetic  vector  potentials, 


A*(z  >i,)  =  i  X  =  i  X  <5'5a) 

p,q  =  -oa  p,q  =  -oo 

oo  oo 

F‘  b  >  z,)  =  i  X  ^  L  ‘<lWp‘l  ■  <5.5b) 

p.q  =  -«  p,q  = 

where  and  /?™  are  complex-valued  scalar  constants.  An  expression  for  the  electric 
field  is  obtained  by  differentiating  Eq.  (5.5)  according  to. 


=  -V  X  F"  -jco^A"  +  V(V-A")  , 

^  JCOSq 


or,  expanding  this  in  Cartesian  coordinates. 


p-<i 

?■<! 

^  =  I 


_  nTE  _  j^pq^pq  pTM 
JPpq^pq  ^pq 


pq 


,  nTE  j^pq^pq  t>TM 

^J^<Spq  -  -f^^pq 


pq 


p-q 


I 


'pq 


Cy^pq- 


(5.6a) 


(5.6b) 


(5.6c) 


and  are  known  quantities,  given  by  Eqs.  (5.1)  and  (5.4).  Equations  (5.6a)-(5.6b) 
are  therefore  a  system  of  two  equations  in  the  two  unknowns  and  Solving  for 

the  unknowns,  the  result  is 
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TE  _  j^Ppq^^pqx  ^pq^^pqy  ^OOy^Qp^Qq^^ 

TM  _  j^^q^^pqx  ■*■  ^OOx^Op^q^  "*■  Ppq^^pqy  "*■  ^00>^0p*5o^)] 


(5.7a) 


(5.7b) 


The  reflection  coefficients  are  defined  as  /?^for  TE  incidence  and  Rqq  for  TM  incidence. 

For  a  strip  grating,  the  reflection  coefficients  are  defined  in  the  same  way.  The 
electric  field  is  first  expanded  as  a  discrete  spectrum  of  plane  waves  of  the  form, 

ee 

E'(z  >  z,)  =X  (e„  +  "  V> '  r.-,  (5.8) 


where 


2mn 

«m  =  . 


(5.9a) 


(5.9b) 


At  z  =  Zf,  the  tangential  electric  fields  due  to  the  induced  current  J  are  represented 
symbolically  as 


Ex-i2  =  z,)  -4(2  =  0) 


(5.10) 


Ey-iZ  =  Zt)  J  iGyx  G^J,  J[7y(z  =  0) 


where 


j(a,)3)  =  ■^Jo( a)  ^  5{a  -  aj  dip  -  kyi) .  (5.11) 


Matching  the  two  representations  at  z  =  z,. 


X  Sia  -  aJSip  -  kyP 


47r  "V  - 

=  —  2j  Gia^.ky,)jQia^)5ia-aJ5ip-kyi),  (5.12) 
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from  which  it  follows  that 

- 

~  ^  ^  •  (5.13) 


If  the  scattered  field  is  assumed  to  be  composed  of  TE  and  TM-polarized  plane- wave  fields 


derived  from  A'^  and  F'^,  where 


A^(z 

'V  dTM.  r 

=  ^Z^Rm  ¥m 

(5.14a) 

m  =  -<» 

m  =  -oo 

F"  (z 

oo 

oe 

(5.14b) 

fn  ~  — 

m  =  -o® 

then 


TE  y[^yi(^ffu  ^Ox^Qrr^  ^rrS^my  ^Qy^Qm)^ 


Rm  = 


ici+ic^i) 


D  ^  — 


(5.15a) 


(5.15b) 


As  before,  and  are  defined  to  be  the  reflection  coefficients  for  TE  and  TM 


incidence,  respectively. 

If  the  incident  field  is  derived  from  unit-amplitude  vector  potentials,  the  reflection 
coefficients  defined  here  have  the  following  properties; 

1.  The  magnitudes  of  Rp^  and  Rp^^  may  be  greater  than  one,  panicularly  for 
fields  scattered  from  the  TM  polarization  into  the  TE  polarization. 


2.  However,  conservation  of  power  requires  that 

I  1  <  1,  for  TE  incidence,  and 
1  /?™  I  <  1,  for  TM  incidence. 


3.  The  reflection  coefficients  for  a  free-standing  perfectly-conducting  plane  at  z  =  0 
are 

R^  =  - 1 ,  for  TE  incidence,  and 

TM 

Rqq  =  +1,  for  TM  incidence. 


67 


5.2  Reflection  Coefficient  for  Finite  Surfaces 

The  reflection  coefficients  defined  for  periodic  surfaces  were  carefully  chosen  to  be 
compatible  with  the  definitions  in  this  section.  In  the  limit  as  the  number  of  elements  in  a 
finite  frequency  selective  surface  approaches  infinity,  the  reflection  coefficient  defined  for 
the  finite  surface  must  be  equivalent  to  the  reflection  coefficient  defined  for  the  periodic 
surface  in  order  to  make  a  direct  comparison  between  the  two. 

Rather  than  attempt  to  justify  each  individual  assumption,  it  is  simpler  to  postulate  a 
reasonable  definition  for  the  reflection  coefficient  and  verify  numerically  that  it  has  the 
desired  properties.  Since  the  reflection  coefficient  for  a  perfectly-conducting  plane  is  equal 
to  ±1,  then,  except  for  a  minus  sign  for  the  TE  polarization,  the  reflection  coefficient 
derived  in  the  previous  section  can  be  interpreted  as  the  complex  amplitude  of  the  plane- 
wave  field  scattered  by  the  surface  in  the  specular  direction  divided  by  the  complex 
amplitude  of  the  plane-wave  field  reflected  from  a  perfectly-conducting  plane.  This  is 
essentially  the  definition  that  will  be  used  for  the  finite  array,  with  the  language  modified  to 
account  for  the  finite  nature  of  the  scatterer. 

The  scattered  field  can  be  viewed  as  a  superposition  of  two  parts:  E^,  the  field 
scattered  from  the  surface  of  the  dielectric  medium  with  the  printed  surface  removed,  and 
E^,  the  field  due  to  the  conduction  current  on  the  printed  surface,  radiating  in  the  presence 
of  the  dielectric  medium.  With  the  printed  surface  removed,  the  plane-wave  field  reflected 
from  the  infinite  dielectric  stratification  E^  is  of  the  form 

E5(z  2  z,)  =  (5,16) 

where  the  vector  constant  is  computed  by  the  transmission-line  method  using  the 
techniques  discussed  in  Chapter  3.  For  a  finite-width  frequency  selective  surface,  this  part 
of  the  scattered  field  is  approximately  proportional  to  the  field  produced  by 


:J,(^,:>..z,)  =  1  e'oo, 


f  -^eT. 

NeTA 

III 

1  2 

2  J 

y\  2  ’ 

2  J 

(5.17) 
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which  is  a  physical  optics  approximation  to  the  tangential  fields  at  the  surface  of  the 
dielectric.  Ng  is  the  number  of  elements  in  a  unit  row  across  the  finite  dimension  of  the 
array;  and  Ty  are  the  x  and  y  periodicities  of  the  printed  surface  (T^  =  and  Ty  = 
T^sinQ.).  The  truncation  factors  and  Py(a,b)  are  zero  except  on  [a,b]  where  they  are 
equal  to  one.  The  far  field  radiated  by  this  surface  field  is  computed  using  standard 
techniques. 

The  remaining  component  of  the  scattered  fields  is  due  to  the  conduction  current  on 
the  printed  surface.  To  keep  the  proper  proportionality  with  the  first  term,  is  the  field 
due  to  a  unit  row  of  elements  radiating  in  the  presence  of  the  dielectric  support  and  is  given 
by  a  spectral-domain  Green's  function  of  the  form. 


11 

_ 1 

I - 

H 

- 1 
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o 

_ i 

Ey.(2  =  Zi)  _ 

.^yx 

^yy. 

y/2  =  o)_ 

(5.18) 


The  far  fields  produced  by  this  equivalent  magnetic  current  sheet  are  obtained  by  evaluating 
the  Fourier  transforms  of  the  electric  field  at  arguments  which  are  determined  by  the 
observation  angles  i9,(p).  If  the  incident  plane-wave  field  is  incident  from  the  direction 
(6i,(pi),  the  reflection  coefficient  is  defined  as  the  total  far- zone  scattered  field  +  E^  in 
the  specular  direction,  {-di.cpi),  normalized  by  E^,  the  far-zone  physical  optics  field 
scattered  by  a  perfectly-conducting  plate  of  area  x  Ty. 


„  £5,(-e„ 0i)  +  O-e.. 

^  "j'g  incidence 


<pi) 

£^e(-g„0,)-hE^g(-g„0.) 

<Pi) 


,  TM  incidence. 


(5.19a) 


(5.19b) 


It  is  easily  verified  that,  for  TE  incidence,  the  physical  optics  approximation  for  the  current 
on  a  perfectly-conducting  plate  produces  a  far-zone  electric  field  that  is  entirely  in  the  (|) 
direction,  and  for  TM  incidence,  the  far  field  is  entirely  in  the  ^  direction. 
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This  definition  is  easily  extended  to  surfaces  which  are  finite  in  both  dimensions,  and 
the  same  definition  can  also  be  used  for  finite  strip  gratings,  with  the  slight  inconvenience 
of  having  to  define  an  arbitrary  period  “7^.”  In  addition,  the  following  properties  have 
been  verified  numerically:  if  the  induced  current  satisfies  the  periodic  condition,  i.e.,  if, 
instead  of  the  current  on  a  finite  frequency  selective  surface,  the  current  on  a  number  of 
elements  of  a  periodic  surface  is  substituted  for  Jo(x,y),  then,  (1)  the  reflection  coefficient 
is  independent  of  the  number  of  elements  A/^,  and  (2),  except  for  a  minus  sign  for  the  TE 
polarization,  the  value  of  the  reflection  coefficient  is  identical  (to  at  least  four  significant 
figures  in  all  our  sample  calculations)  to  the  value  obtained  by  using  the  plane-wave 
definition. 


5.3  Finite  Strip  Grating 

In  Chapter  2,  it  was  shown  that  for  a  finite  planar  strip  grating  illuminated  by  an 
arbitrary  plane  wave,  the  current  density  J  on  the  infinite  plane  could  be  represented  in 
terms  of  the  current  density  Jo(a:)  along  the  cut  y  =  0.  In  the  transform  domain, 

)  =  2k  Jo(a)5()3  -  ky) . 

In  Chapter  3,  an  integral  equation  for  the  fundamental  unknown  Jo(x  )  was  derived,  and 
the  spectral-domain  method  of  moments  was  applied  to  generate  a  matrix  equation;  i.e.,  if 

/ 

J=  1 


then 


*00  00 

f  j'  G{a,ky)  j/a)e^Vrfa  =  f  E‘"^(a:,y)  tic ,  i  =  1,  2,  ...,  /  . 

A  =-00 


In  Chapter  4,  the  essential  theoretical  and  numerical  considerations  were  discussed  that 
pertain  to  implementing  the  spectral-Galerkin  equation  efficiently  and  accurately  on  the 
computer.  The  only  remaining  steps  involve  choosing  a  linearly-independent  basis  {jy}  for 
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expanding  the  unknown  and  solving  the  resultant  system  of  equations  for  the  coefficients 

{Cj). 

The  transverse  current  7^  is  discretized  using  the  trigonometric  series 
r  kTT  1 

=  sin  (Ct  -  Xi)  +  h)  \Px(xk  -  h,X;^-h  h),  k  =  \,  2,  .  .  (5.20) 

and  the  parallel  component  of  the  current  JyQ  is  represented  by  functions  of  the  form, 

jykM  - - ,  .  —PxiXk  -h,Xk  +  h),  k=  2,  ,  (5.21) 


ix-Xk) 


where  7'n(^)  is  the  Chebychev  polynomial  of  order  n, 

r„(jr)  =  cosin  c^  '  c)) , 

h  is  the  half-width  of  the  strips  (h  =  a/2),  x/^  is  the  position  of  the  basis  function  relative  to 
the  origin,  and  the  truncation  factor  P^,  which  limits  the  extent  of  the  functions  to  the 
surfaces  of  the  conductors,  is  defined  as 


(  0,  otherwise 

Note  that  each  term  in  the  expansion  of  JyQ  has  a  built-in  edge  singularitv  to  model  the 
physical  behavior  of  the  current.  The  first  few  terms  of  this  series  are  shown  in  Fig.  5.2. 
The  Fourier  transforms  of  the  basis  functions  are  [34] 


jxkioc)  =  -jh  2  )|  sinc^-^  -  -  (-l)^sinc^-^  +  a/^e  (5.22) 

jyi^ia)  =  i-}f~^ nh  /(i_i)(-a/j)  ,  (5.23) 

where  y  =  V-1,  sinc(jc)  =  sin(a:)/x,  and  J^ix)  is  the  nth-order  Bessel  function  of  the  first 
kind. 
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Figure  5.2.  Chebychev  polynomials  modified  to  incorporate  the  edge 
singularity. 


To  verify  that  these  basis  and  testing  functions  have  sufficient  smoothness  to 
guarantee  convergence  of  the  inner  products,  the  asymptotic  behavior  of  the  product 
}r  icc)Gia,kyi)ij(a)  is  examined  for  each  combination  of  basis  and  testing  function,  i.e., 

^«(«)  =  Jxk(odG^^ia,ky)j^k(cc) 

Rxyioc)  =  hk^a)G^y{a,ky)jyf,{a) 

Ryyia)  =  jyk(a)Gyyia,ky)  jyi^ia)  . 

The  two  sin(a)/a  terms  combine  in  such  a  way  that  the  net  convergence  of  jxk((^)  is  not 

1/a,  but  1/a^.  This  might  have  been  expected,  since  the  trigonometric  basis  functions  have 

the  same  degree  of  differentiability  as  the  triangular  basis  function,  whose  Fourier 

transform  is  which  is  well  known.  The  asymptotic  expression  for  J„(x)  is  given 

a 

by 
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If  each  term  in  the  product  jj'*"(a)G( a,/:y,)jy( a)  is  replaced  by  its  asymptotic  form,  then, 
except  for  a  constant,  the  asymptotic  expressions  for  the  integrands  R^,  Rjy,  and  Ryy  are 
given  by 


Rj^icx)  — > 
Rj^yicc)  — > 

Ryyicc))  - > 


(5.24a) 

(5.24b) 

(5.24c) 


For  this  choice  of  basis  and  testing  functions,  the  integrands  and  summands  converge  at 
least  as  fast  as  i/oP'. 


The  inner  products  formed  by  applying  the  spectral-domain  method  of  moments 

contain  either  infinite  summations  or  integrals  with  infinite  limits  of  integration.  If  no 

extrapolation  techniques  are  used  to  predict  the  contribution  of  the  tail  regions  [35],  the 

summations  or  integrations  are  truncated  to  a  finite  area  of  the  a-/3  plane.  When  truncating 

the  spectrum,  however,  the  limits  of  summation  or  integration  must  encircle  the  highest 

spectral  component  of  the  current;  otherwise,  the  moment  matrix  will  be  ill-conditioned 

[36].  For  the  trigonometric  series,  the  peak  of  the  highest  spectral  component  of  the 

current  occurs  at  a  =  .  Therefore,  the  first  opportunity  to  truncate  the  spectrum  is  at 

TC 

the  first  zero  beyond  the  peak  or  at  {Nsin  +  2)^.  By  a  similar  argument,  the  spectral 

K 

bandwidth  of  the  Chebychev  series  is  approximately  {Ncheb  +  i)^-  in  all  of  the  results  in 
this  section,  the  spectrum  has  been  truncated  at  twice  the  maximum  of  these  two  numbers. 

Electromagnetic  scattering  from  a  free-standing  periodic  strip  grating  has  been  solved 
exactly  by  Weinstein  [10]  using  the  generalized  Wiener-Hopf  technique.  Figure  5.3  gives 
an  idea  of  the  accuracy  of  the  numerical  results  by  comparing  the  magnitude  of  the 
reflection  coefficient  computed  by  Weinstein’s  method  with  the  results  of  applying  this 
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theory  to  the  same  problem.  The  incident  electric  field,  which  is  incident  normally  on  the 
surface,  is  polarized  either  parallel  or  perpendicular  to  the  axes  of  strips.  The  numerical 
results  have  been  computed  assuming  four  trigonometric  basis  functions  and  four  modified 
Chebychev  polynomials  to  represent  the  unknown.  Since  the  incident  angles  {di,(pi)  = 
(0,0)  are  not  allowed  when  computing  the  incident  field  by  the  vector  potential  method, 
{01, (pi)  =  (0.1°,  0.1°)  have  been  used  for  the  numerical  calculations.  The  worst-case 
deviation  of  the  results  is  4.3  percent. 

To  illustrate  the  effect  of  truncation.  Figs.  5.4  and  5.5  show  the  frequency  response 
of  the  periodic  strip  grating  of  Fig.  5.3,  truncated  to  3,  7,  and  15  strips.  This  case  was 
considered  by  Cwik  and  Mittra  [9]  using  the  equivalent  of  64  subsectional  basis  functions 
per  strip,  for  the  two  polarizations  of  the  incident  field.  In  the  15-strip  example,  this 
requires  filling  a  960  x  960  matrix.  The  results  of  Figs.  5.3  and  5.4,  which  are  in  excellent 
agreement  with  the  results  of  Cwik  and  Mittra,  were  computed  using  8  basis  functions  per 
strip,  or  by  filling  a  120  x  120  matrix.  This  underlines  the  importance  of  choosing  a  basis 
set  that  is  capable  of  modeling  the  current  in  a  minimum  number  of  terms. 

As  pointed  out  by  the  previous  investigators  [9],  the  effect  of  finite  dimensions  is 
greater  when  the  incident  electric  field  is  polarized  parallel  to  the  truncation  boundary,  i.e., 
when  the  dominant  component  of  the  induced  current  is  flowing  parallel  to  the  edge  of  the 
finite  array.  In  Fig.  5.4,  where  the  incident  electric  field  is  polarized  parallel  to  the 
truncation  boundary,  the  fields  scattered  by  a  finite  strip  grating  can  be  significantly 
different  from  the  fields  produced  by  an  equal  number  of  strips  in  an  infinite  array;  this 
difference  becomes  less  as  the  number  of  strips  is  increased.  Conversely,  when  the 
incident  electric  field  is  polarized  perpendicularly  to  the  edge  of  the  array  over  this 
frequency  range,  the  difference  between  the  periodic  solution  for  the  reflection  coefficient 
and  the  finite  solution  is  minimal,  even  for  as  few  as  3  strips.  Because  of  the  properties  of 
the  reflection  coefficients  mentioned  in  Sections  5.2  and  5.2.  the  close  agreement  in  Fig. 
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Figure  5.3.  Comparison  of  numerical  and  exact  reflection  coefficients  for  a 
periodic  strip  grating:  strip  width  a  =  0.5  cm,  periodicity  b  = 
1  cm.  (a)  polarized  parallel  to  the  axis  of  the  strips, 
(b)  polarized  perpendicular  to  the  axis  of  the  strips. 
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Figure  5,4  (a)  Magnitude  and  (b)  phase  of  the  computed  reflection 

coefficient  for  a  plane  wave  incident  on  several  free-standing 
strip  gratings  of  3,  7,  and  15  strips,  compared  to  the  periodic 
result:  TE  polarization,  strip  width  =  0.5  cm,  periodicitv  = 
1  cm,  =  (0.1°,  0,1°). 
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Figure  5  5.  (a)  Magnitude  and  (b)  phase  of  the  computed  reflection 
coefficient  for  a  plane  wave  incident  on  several  free-standing 
strip  gratings  of  3,  7,  and  15  strips,  ccmpared  to  the  periodic 
result;  TM  polarization,  stnp  width  =  0.5  cm,  periodicity  = 
1  cm,  =  (0.1°,  0.1°). 


77 


I 

e 


5.5  indicates  that,  for  perpendicular  polarization,  the  field  scattered  in  the  specular  direction 
by  a  finite  number  of  strips  can  be  accurately  approximated  by  the  field  produced  by  an 
equal  number  of  strips  whose  currents  have  been  obtained  from  a  periodic  model. 

One  of  the  reasons  why  the  periodic  approximation  is  valid  for  a  finite  array  is  that  the 
edge  effect  is  localized.  This  is  illustrated  in  Figs.  5.6(a)  and  5.6(b).  A  periodic  surface  is 
characterized  bv  the  fact  that  the  induced  current  on  any  one  cell  differs  from  the  current  on 
any  ot.her  cell  only  by  the  uniform  progressive  phase  shift  of  the  incident  field.  Therefore, 
the  magnitudes  of  the  current  densities  on  a  number  of  cells  will  be  identical  in  both  shape 
and  amplitude.  If  the  finite-array  problem  is  solved  using  the  same  approximations  that 
were  used  in  the  periodic  problem,  then  any  deviation  in  the  behavior  of  the  current  from 
this  uniform  amplitude  and  shape  is  entirely  due  to  the  edge  effect.  A  plot  of  the  magnitude 
of  the  current  density  across  a  number  of  cells  of  a  15-element  strip  grating  is  shown  in 
Fig,  5.6(a)  for  parallel  polarization  and  Fig.  5.6(b)  for  perpendicular  polarization.  Only 
half  of  the  cell  currents  are  plotted  since  the  incident  field  and  geometry  are  symmetrical. 
Also  shown  are  the  corresponding  currents  on  a  periodic  strip  grating.  From  these  results, 
it  is  evident  that  there  is  a  locally-confined  edge  region  that  penetrates  only  about  three 
strips  into  the  interior  of  the  array.  As  the  array  becomes  larger,  the  edge  region  occupies  a 
smaller  proportion  of  the  total  area.  The  interior  region,  however,  where  the  currents 
closely  match  those  of  the  periodic  array,  occupies  a  larger  proportion  of  the  total  area. 

5.4  Crossed  Dipole  Element 

The  crossed  dipole  element  was  originally  intended  to  be  a  simple  extension  of  the 


linear  dipole.  It  had  been  known  for  some  time  that  for  an  incident  plane-wave  field  of  the 
proper  polarization  and  frequency,  a  periodic  array  of  thin  conducting  dipoles  was  capable 
of  reflecting  the  field  as  though  it  were  a  sofd  sheet  of  perfectly  conducting  metal  [37]. 
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Figure  5.6.  Dominant  currents  induced  on  a  15-element  finite  strip  grating 
along  the  cut  >’  =  0:  (a)  parallel  polarization,  (b)  perpendicular 
polarization.  Strip  width  =  0.5  cm,  periodicity  =  1  cm, 
=  0.1°). 
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The  crossed  dipole  was  an  attempt  to  make  this  effect  independent  of  the  polarization  of  the 
incident  field. 

The  crossed  dipole  has  been  a  popular  choice  in  applications  requiring  circular 
polarization.  For  certain  skew  angles,  the  crossed  dipole  array  has  the  property  that  at 
normal  incidence,  the  surface  looks  the  same  to  each  of  the  two  linear  components  of  a 
circularly-polarized  incident  field.  On  reflection  or  transmission,  the  relative  amplitudes 
and  phases  of  the  two  components  are  preserved,  and  the  scattered  fields  are  circularly- 
polarized  as  well.  Examples  of  other  elements  having  the  same  property  are  Jerusalem 
crosses,  square  patches,  and  rings.  But  because  of  the  ease  of  fabricating  crossed  dipole 
arrays  and  the  apparent  simplicity  of  modeling  them  analytically,  most  of  the  experimental 
and  theoretical  work  has  been  conducted  with  crossed  dipole  arrays  [8,  14,  17,  38-40]. 

However,  the  apparent  simplicity  of  modeling  periodic  arrays  of  crossed  dipoles  was 
deceptive.  This  was  due  to  the  fact  that  the  current  induced  on  the  two  perpendicular  arms 
of  the  cross  was  not  simply  the  superposition  of  the  currents  on  two  isolated  linear  dipoles. 
In  addition  to  these,  another  term  was  present,  sometimes  referred  to  as  the  “crooked 
mode,”  due  to  the  mutual  coupling  between  the  two  arms.  In  1977,  King  published  his 
solution  for  the  zero-order  currents  on  the  cross  [41],  which  included  the  new  term,  and 
soon  afterward,  Pelton  and  Munk  showed  that  it  was  the  presence  of  this  term  that  caused 
the  anomalous  antiresonance  that  had  been  observed  in  the  frequency  response  of  crossed 
dipole  arrays  for  oblique  incidence  [39]. 

The  problem  was  later  revisited  by  Tsao  and  Mittra  [8],  who  incorporated  the  so- 
called  “crooked  mode”  by  using  a  special  junction  basis  function.  The  basis  functions  used 
here  are  the  same  as  those  used  by  Tsao  and  Mittra  in  their  treatment.  Referring  to  Fig. 
5.7,  if  L  is  the  length  and  W  is  the  width  of  each  of  the  two  dipole  members,  then  the  basis 
functions  can  be  expressed  in  the  form 

.  ,  ,  -IkKf  L\\„(  -L  L\^(  -W  W'\  , 

AAy)  =  sin|^— J  .  /t=  1,  2,  .  .  (5.25a) 
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jykf<^.y)  =  ^  2,  .  .  .,  (5.25b) 

jjix,y)  =  X  sgn(A:)  cos(^^)p^(^-Y,Yy3,(^-|^, y) 

-  y  sgn(y)  j  j  ,  (5.25c) 

where  only  the  basis  functions  at  the  origin  have  been  shown.  The  truncation  factors 
Px'a.h)  and  Py‘[a,b)  are  defmea  to  oe  zero  except  on  the  interval  [a,b],  where  they  are  equal 
to  one. 


Figure  5.7.  Geometry  of  the  crossed-dipole  array. 


The  Fourier  transforms  of  the  basis  functions  are 
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(5.26c) 


The  junction  function  is  only  included  to  provide  a  minimum-term  '•epresentation  for 
the  current  on  the  cross.  Basis  functions  of  either  subdomain  or  entire-domain  type  are 
capable  of  modeling  the  current  at  a  junction,  provided  the  number  of  terms  is  sufficiently 
high.  By  using  the  junction  function,  the  object  is  to  reduce  the  number  of  unknowns  from 
approximately  10  per  wavelength  to  approximately  2.5  per  wavelength. 

These  basis  functions  give  the  weakest  form  of  convergence.  Asymptotically,  the 
Fourier  transform  of  the  junction  basis  function  decays  as  ^  ^ .  Combining  this  with  the 
asymptotic  expression  for  the  spectral  Green's  function,  the  product  j/’i’Gjy  decays  only  as 
fast  ^  ^  for  the  junction-junction  interaction.  If  p  is  held  constant  while  an  integration 
is  performed  over  a,  which  is  the  procedure  suggested  by  Eq.  OAl),  the  inner  product 
will  not  converge.  However,  in  polar  coordinates  [42],  where  a  =  kpcoscp  and  P  =  kpSincp, 
the  worst-case  convergence  rate  is  1/k^sin  <p,  i.e.,  except  for  along  the  radius  (p  =  0,  the  kp 
integration  converges. 

Sample  results  using  this  basis  set  are  shown  in  Fig.  5.8,  for  a  geometry  similar  to 
one  in  the  literature  [15].  The  geometry  is  a  periodic  array  of  cross-shaped  conducting 
patches  (T-rj^  =Tri2  =  ^  cm,  L  =  0.6875  cm,  W  =  0.0625  cm,  Q  =  90  deg)  printed  on  a  3 

mm  dielectric  substrate  of  relative  permittivity  e^.  =  4  -  j.02.  The  frequency  response  has 
been  computed  twice,  once  using  rooftop  basis  functions  (12  per  arm),  and  once  using  the 
entire-domain  functions  of  Tsao  and  Mittra  (2  per  arm  plus  the  junction  function).  The 
agreement  between  the  results  is  quite  good.  It  is  interesting  to  note  that  the  addition  of  a 
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small  amount  of  loss  in  the  dielectric  region  prevents  the  reflection  characteristic  from 
reaching  a  maximum  of  1 .0  at  resonance. 

If  the  same  element  and  lattice  are  assumed  for  frequency  selective  surfaces  of  3,  7, 
and  15  patches  wide,  most  of  the  difference  in  the  plot  of  the  reflection  coefficient  occurs 
near  the  resonant  frequency.  Figures  5.9  and  5.10  show  the  variation  in  the  reflection 
coefficient  for  different  array  sizes.  As  before,  the  greatest  difference  between  the  finite 
and  periodic  results  occurs  for  the  smallest  number  of  elements,  the  difference  becoming 
less  as  more  elements  are  considered.  When  the  dominant  component  of  the  current  is 
flowing  parallel  to  the  edge  of  the  array  (TE  polarization),  the  edge  has  more  influence  on 
the  reflected  fields;  when  the  dominant  current  is  flowing  perpendicularly  to  the  edge  of  the 
array,  the  edge  effect  is  less  significant. 

As  a  second  example,  the  free-standing  surface  of  Pelton  and  Munk  [39]  is 
considered,  which  is  periodic  in  skewed  coordinates:  T’tjj  =  2.44  cm,  7'^2  =  i-725  cm,  L 
=  1.51  cm,  W  =  0.044  cm,  =  45  deg.  Figure  5.11  compares  computed  values  of  the 
reflection  coefficient  from  a  periodic  surface  and  their  measured  results  for  oblique 
incidence  (0/  =  30°),  showing  the  characteristic  antiresonance  of  the  crossed-dipole  array. 
For  this  example,  two  trigonometric  basis  functions  per  arm  plus  the  junction  function  were 
used  to  approximate  the  current  on  the  cross.  The  Floquet  harmonics  were  summed  out  to 
a  point  which  was  twice  the  spectral  bandwidth  of  the  highest  (trigonometric)  component  of 
the  current.  Four  basis  functions  per  arm  reproduced  their  results  exactly. 

If  the  periodic  surface  of  Fig.  5. 1 1  is  truncated  to  15  elements  wide,  the  effect  on  the 
reflection  coefficient  is  shown  in  Fig  5.12.  The  reflection  coefficient  for  the  periodic 
surface  is  shown  for  comparison.  Both  results  were  computed  using  the  same 
approximations:  two  trigonometric-type  basis  functions  per  arm  plus  the  junction  function, 
encircling  that  portion  of  the  spectrum  that  was  less  than  twice  the  spectral  radius  of  the 
(trigonometric)  basis  functions. 
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Figure  5.8.  (a)  Magnitude  and  (b)  phase  of  the  reflection  coefficient  for  a 
periodic  array  of  cross-shaped  elements  on  a  dielectric 
substrate,  computed  with  two  different  sets  of  basis  functions. 
=  T’tjo  =  1  cm,  L  =  0.6875  cm,  W  =  0.0625  cm,  n  = 
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Figure  5.9.  Computed  values  for  the  reflection  coefficient  for  finite-width 
frequency  selective  surfaces  of  3,  7,  and  15  patches  on  a 
dielectric  backing.  Same  element  and  lattice  as  in  Fig.  5.8, 
TE  polarization  =  (0.1°,  0,1°): 

(a)  frequency  range  8-18  GHz,  (b)  12-14  GHz. 
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Figure  5.10.  Computed  v^’oes  for  the  reflection  coefficient  for  finite-width 
frequency  selective  surfaces  of  3,  7,  and  15  patches  on  a 
dielectric  backing.  Same  element  and  lattice  as  in  Fig.  5.8, 
TM  polarization  (E'"*^  =  =  (0.1°,  0.1°); 

(a)  frequency  range  8-18  GHz,  (b)  12-14  GHz. 
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Figure  5.11.  Comparison  of  computed  and  measured  data  for  the 
frequency  selective  surface  of  Pelton  and  Munk  [39]; 
periodic  array,  TM  polarization  {di,(pi)  = 

(30^  0°). 
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Figure  5.12.  Computed  reflection  coefficients  for  a  15-element  frequency 
selective  surface  of  finite  width  and  its  periodic  counterpart; 
TM  polarization  (E*"c  (di,(pi)  =  (30°,  0°). 
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These  results  show  the  error  in  using  the  periodic  approximation  to  predict 
electromagnetic  scattering  from  a  finite  array  of  various  sizes,  when  the  frequency  of  the 
incident  field  is  near  the  first  resonance  of  the  element.  If  a  finite-width  array  is  15  patches 
wide,  the  error  in  the  scattered  field  using  the  periodic  approximation  is  on  the  order  of  two 
per  cent.  However,  if  the  frequency  is  much  lower  than  the  first  resonant  frequency  of  the 
element,  a  low-frequency  resonance  can  be  excited  on  the  array,  as  shown  in  Figs.  5.13 
and  5.14,  even  when  there  is  no  direct  electrical  connection  between  the  two  edges.  In 
both  of  these  examples,  the  array  is  embedded  between  two  identical  dielectric  layers  of 
relative  permittivity  4  -  j.Q2  and  thickness  r  =  .5  cm;  In  the  first  example,  the  array  is 
exactly  two  wavelengths  wide  at  the  frequency  of  the  incident  field.  The  two-wavelength 
envelope  is  plainly  visible  in  the  magnitude  plot.  In  the  second  example,  the  array  is  four 
wavelengths  wide  at  the  frequency  of  the  incident  field.  Since  the  edge  effect  is  not 
localized  in  this  frequency  range,  approximating  the  fields  from  a  periodic-array 
formulation  would  be  significantly  in  error.  Such  models  cannot  predict  the  interaction  of 
two  edges,  which  is  clearly  present  in  Figs.  5.13  and  5.i4. 

None  of  the  previous  examples  illustrate  the  versatility  of  the  multilayer  Green’s 
function  developed  in  Sec.  3.1.  As  a  final  example,  sample  results  for  a  periodic  surface 
embedded  in  a  four-layer  dielectric  medium  are  presented.  The  geometry  is  shown  in  Fig. 
5.15,  which  was  described  by  Contu  and  Tascone  [43].  The  element  is  a  crossed  dipole, 
whose  lattice  and  dimensions  are  given  by  Tiji  -  0.  65  cm,  Trj2  =  0-46  cm,  L  =  0.555  cm, 
IT  =  0.015  cm,  42  =  45  deg,  as  defined  in  Fig.  5.7.  The  results  are  presented  in  terms  of 
the  magnitude  of  the  reflection  coefficient,  which,  for  frequencies  low  enough  that  only  one 
propagating  mode  exists  in  the  free-space  regions  above  and  below  the  surf..ce,  can  be 
accurately  approximated  by 
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Figure  5. 13  Low-frequency  resonance  on  a  finite  array  of  15  cross-shaped 
conducting  patches  embedded  between  two  identical  dielectric 
layers.  Tt^i  =Tri2=  1-6  cm,  L  =  1.51  cm,  W  =  0.044  cm,  Q 
=  90  deg,  £r,,  C;.,  =  4  -  j.02,  t\,  t2  =  0.5  cm,  frequency  = 
2.5  GHz,  (6i,(pi)  =  (1°,  1°);  (a)  magnitude,  (b)  phase. 
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Figure  5.14.  Low-frequency  resonance  on  a  finite  array  of  15  cross-shaped 
conducting  patches  embedded  between  two  identical  dielectric 
layers.  Ttj^  =  7’tj2  =  cm,  L  =  1.51  cm,  W  =  0.044  cm,  Q 
=  90  deg,  £r,,  £r2  =  ^  -  J-02,  ti,  t2  =  0.5  cm,  frequency  = 
5  GHz,  idi,(pi)  =  (1°,  1°):  (a)  magnitude,  (b)  phase. 
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Figure  5.15.  Four-layer  dielectric  medium  (Contu  and  Tascone  [43]). 


The  numerical  results  of  Fig.  5.16  were  computed  using  four  trigonometric-type 
basis  functions  on  each  of  the  two  anns  plus  the  junction  function.  For  the  interactions  that 
do  not  involve  the  junction  basis  function,  the  inner  products  were  evaluated  by  truncating 
the  spectra  at  eight  times  the  highest  spectral  bandwidth  of  the  trigonometric  basis 
functions.  For  those  matrix  elements  that  involve  the  junction  basis  function,  the  inner 
products  were  evaluated  by  summing  out  to  sixteen  times  the  highest  spectral  bandwidth  of 
the  trigonometric-type  basis  functions. 

The  discrepancy  between  the  two  results  at  higher  frequencies  can  be  reduced  for 
the  TE  polarization  by  including  more  basis  functions  in  the  expansion  of  the  unknown. 
However,  this  does  not  explain  the  discrepancy  for  the  TM  polarization,  for  which  the 
addition  of  up  to  twenty  basis  functions  per  arm  does  not  improve  the  agreement  between 
the  two  results.  The  accuracy  of  these  results  has  also  been  confirmed  by  an  independent 
program  using  twenty  rooftop  basis  functions  per  wavelength.  The  difference  might  be 
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Figure  5.16.  Computed  values  of  the  transmission  coefficients  for  the 
geometry  of  Contu  and  Tascone  [43]:  {6,(p)  =  (30°,  45°). 
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explained  by  one  of  the  following  facts.  (1)  Unlike  the  paper  by  Contu  and  Tascone,  no 
attempt  was  made  here  to  incorporate  the  edge  singularity  in  the  representation  for  the 
current.  This  is  usually  not  a  significant  source  of  error  for  narrow  elements.  (2)  Since 
Contu  and  Tascone  have  employed  the  scattering-matrix  approach  to  incorporate  the 
dielectric  support,  perhaps  the  difference  can  be  blamed  on  the  fact  that  typically  a  large 
number  of  Floquet  space  harmonics  must  be  included  in  the  scattering  matnces  as  the  layers 
become  extremely  thin,  as  they  are  in  this  example. 
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CHAPTER  6 

ARBITRARY  INCIDENCE  AND  NONUNIFORM  ARRAYS 

When  a  frequency  selective  surface  is  considered  to  be  finite,  the  vector-Floquet 
nature  of  the  fields  vanishes,  i.e.,  the  scattered  field  is  no  longer  a  discrete  spectrum  of 
plane  waves,  regardless  of  the  nature  of  the  incident  field.  In  this  case,  it  is  not  necessary 
to  restrict  one's  attention  to  plane-wave  sources,  or  to  truncated  versions  of  periodic 
surfaces.  Arbitrary  incident  fields  can  easily  be  treated  in  the  context  of  finite  arrays,  and 
all  of  the  tools  that  were  developed  to  predict  the  currents  on  uniform  finite  arrays  are 
directly  applicable  to  nonuniform  arrays.  The  designer  is  therefore  capable  of  modeling  the 
effects  of  varying  the  size  and  shape  of  the  elements,  as  well  as  the  array  lattice,  allowing 
more  flexibility  to  achieve  his  design  goals.  The  reasons  for  doing  so  will  become  apparent 
as  we  proceed. 

For  most  of  the  unit  cell  geometries  used  in  practice,  the  null  frequencies  in  the 
reflection  or  transmission  characteristics  are  very  sensitive  to  the  incident  angles  of  the 
incident  field.  Since  FSS  structures  are  typically  narrow  band,  it  is  essential  that  a  surface 
be  designed  for  the  incident  angles  that  exist  in  a  given  system.  However,  in  numerous 
practical  situations,  the  incident  field  is  not  a  plane  wave,  and  an  entire  range  of  incident 
angles  exists  at  the  surface  of  the  array.  This  is  the  case,  for  example,  when  a  frequency 
selective  surface  is  in  the  Fresnel  field  of  an  electromagnetic  horn.  In  this  chapter,  the 
effects  of  a  non-plane-wave  source  will  be  modeled  using  a  very  simple  and  natural 
extension  of  the  plane-wave  situation,  demonstrating  that  a  serious  degradation  in 
performance  results  when  a  realistic  source  is  used.  A  method  is  proposed  to  compensate 
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for  the  effects  of  non-plane-wave  illumination  by  tapering  the  geometry  of  the  surface 
across  the  array. 

The  problem  that  will  be  considered  is  a  dichroic  subreflector  in  an  offset 
configuration.  Because  of  their  ease  of  manufacture  and  analytical  treatment,  planar,  rather 
than  hyperbolic,  surfaces  are  often  used  in  this  application.  In  order  to  isolate  the 
significant  variables  and  to  reduce  the  problem  to  a  manageable  size,  the  analysis  is 
performed  for  a  frequency  selective  surface  of  finite  width,  permitting  the  incident  field  to 
vary  in  a  non-plane-wave  fashion  only  along  the  finite  dimension  of  the  array.  The  x 
variation  of  the  field  is  arbitrary,  but  the  y  variation  is  assumed  to  be  the  phase  dependence 
associated  with  a  uniform  plane  wave.  The  geometry  is  illustrated  in  Fig.  6.1.  In  this 
case,  the  source  is  close  enough  to  the  array  that  the  incident  field  cannot  be  considered  a 
plane  wave. 


Figure  6.1.  A  finite- width,  free-standing  array  of  y-directed  thin  dipoles 
illuminated  by  an  electromagnetic  hom. 


The  application  of  the  spectral-Galerkin  technique  results  in  a  matrix  equation  of  the 


form. 


95 


27tT 


sin^2  1 1  }](a,pj  G («,/?„)  ijia,p^)da 

~J  J  E“'‘'U,}')  dxdy  ,  i=  \,  2,  ...,  1 .  (6.1) 


;  =  1  n  =  ~o 


For  plane-wave  excitarion,  the  kernel  of  the  double  integral  on  the  right-hand  side  of  Eq. 
(6.1)  is  of  Fourier- transform  type,  which  allows  the  integrations  to  be  carried  out 
analytically  (cf.  Sec.  3.5).  This  prompts  us  to  write  an  arbitrary  field  in  terms  of  a  Fourier 
series,  which,  for  the  one-dimensional  case  of  interest  to  us,  takes  the  form 

M 

ErUy)  =  ^  ^  .  (6.2) 

m  =  -M 

where  Ng  is  the  number  of  elements,  is  the  average  period  along  the  x  axis,  and  NgT^  is 
the  width  of  the  array.  The  vector  constants  Cm  are  determined  in  the  usual  manner  by 
multiplying  both  sides  of  Eq.  (6.2)  by  the  complex  conjugate  of  the  exponential  term  and 
integrating  over  the  width,  which  may  be  done  for  all  in  parallel  via  the  fast  Fourier 
transform.  When  this  representation  is  substituted  for  the  incident  field,  the  two- 
dimensional  integrations  on  the  right  can  be  done  analytically,  as  before,  resulting  in  a 
simple  series  expression  for  the  right-hand  side.  It  should  be  emphasized  that  the 
numerically-intensive  process  of  filling  the  moment-method  matrix  is  done  only  once,  and 
the  resultant  system  of  equations  is  solved  using  a  composite  right-hand  side. 

In  order  to  determine  the  effect  of  non-plane- wave  incidence  on  the  scattered  field,  the 
reflection  coefficient  is  defined  as  the  far-zone  scattered  electric  field  in  the  specular 
direction  due  to  a  unit  row  of  dipoles,  normalized  with  respect  to  the  physical  optics 
scattering  of  a  perfectly  conducting  sheet  of  the  same  physical  area.  Since  the  incident  field 
has  been  decomposed  into  a  spectrum  of  plane  waves  according  to  Eq.  (6.2),  the 
calculation  is  easily  performed.  As  in  the  periodic  case,  the  magnitude  of  the  reflection 
coefficient  reaches  a  maximum  of  1.0  at  the  element  resonances. 
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As  an  example,  the  x  variation  of  the  incident  field  is  assumed  to  be  the  H-plane  cut 
of  the  far-field  pattern  of  an  E-plane  sectoral  horn.  It  is  important  that  the  reader  realize  that 
the  far-field  pattern  was  chosen  only  for  convenience  in  order  to  give  a  realistic  model  for 
the  magnitude  and  phase  distribution  of  the  incident  field  on  the  array.  The  horn  model  is 
described  by  Balanis  [44],  in  which  the  magnitude  distribution  of  the  fields  at  the  mouth  of 
the  horn  is  assumed  to  be  the  same  as  the  distribution  for  the  dominant  TEio  mode  in  the 
feed  waveguide,  and  the  phase  distribution  is  assumed  to  be  the  phase  distribution  of  a 
cylindrical  wave  emanating  from  the  apex  of  the  horn  (see  Fig.  6.2).  The  far-field  pattern 
of  this  aperture  field  is  computed  using  standard  techniques,  where  the  integrals  are 
evaluated  in  terms  of  cosine  and  sine  Fresnel  integrals. 


Figure  6.2.  E-plane  horn. 

The  frequency  selective  surface  chosen  for  the  numerical  study  is  a  finite-width 
version  of  the  periodic  surface  described  by  Ott  et  al.  [37].  The  unit  cell  is  defined  in  Fig. 
6.3,  where  Tx  =  Ty  =  1.78  cm,  L  =  1.27  cm,  W  =  0.127  cm,  Q  =  90  deg  Fifteen 
y-directed  thin  dipoles  are  assembled  along  the  x  axis  and  then  replicated  in  a  periodic 
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Figure  6.3.  Unit  cell  in  a  dipole  array. 

fashion  along  the  positive  and  negative  y  axes  to  generate  the  complete  surface.  The 
elements  are  assumed  to  be  thin  enough  that  the  transverse  current  on  each  element  (i.e., 
Jx  )  is  negligible,  as  well  as  the  x  variation  of  Jy.  As  a  result,  the  operator  is  discretized 
using  rooftop  basis  functions  and  Galerkin  testing. 

The  dimensions  of  the  horn  are  a  =  0.5A,  b  =  0.25A,  b\  =  2.75A,  p\  =  6.0A  at 
10  GHz,  which  are  held  fixed  as  the  frequency  is  varied.  The  phase  center  of  the  horn  is 
located  along  the  radial  line  (6,(p)  =  (30°,  180°),  at  a  range  of  26.7  cm  from  the  origin,  as 
defined  in  Fig.  6.1,  and  is  oriented  so  that  the  incident  field  is  parallel  to  the  elements.  For 
these  dimension'',  the  effective  incident  angles  are  0  =  0°  at  the  left  edge,  9  =  30°  at  the 
origin,  and  0  =  49.1°  at  the  extreme  right  edge  of  the  array.  A  plot  of  the  incident  field 
distribution  along  the  x  axis  (magnitude  and  phase)  is  shown  in  Fig.  6.4,  at  a  frequency  of 
10  GHz,  and  assuming  ky  =  0.  The  maximum  field  intensity  occurs  at  a  point  which  is 
midway  between  the  left  edge,  where  the  horn  is  closest  to  the  array,  and  the  origin,  which 
is  along  the  maximum  of  the  horn's  field  pattern.  The  change  in  slope  of  the  phase  plot  is 
due  to  the  change  in  effective  incident  angle  on  the  array. 

The  power  reflection  coefficient  for  a  15-element  array  illuminated  by  this  incident 
field  is  shown  in  Fig.  6.5.  Also  shown  is  the  reflection  curve  for  a  periodic  surface 
iUuminated  by  a  plane-wave  field  incident  from  (0,  ip)  =  (30°,  180°). 


Incident 
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Figure  6.5.  Power  reflection  coefficient  versus  frequency:  plane-wave 
scattering  from  a  periodic  surface,  and  scattering  from  a 
tapered  finite-width  FSS  of  15  elements  due  to  an  incident 
field  that  varies  as  the  H-pIane  cut  of  an  E-plane  sectoral  horn. 


As  is  evident  from  the  figure,  non-plane-wave  illumination  of  the  surface  causes  significant 
distortion  in  its  reflection  curve.  About  ten  per  cent  of  the  available  power  is  lost  at 
resonance  and  about  20  per  cent  of  the  bandwidth. 

The  source  of  the  problem  appears  to  be  the  change  in  effective  illumination  angle  on 
the  array.  Figure  6.6  is  a  plot  of  the  resonant  frequency  versus  incident  angle  for  the 
element  used  in  this  example.  The  resonant  frequency  shifts  from  11.2  GHz  at  0=  0°  to 
9.25  GHz  at  6  =  50°.  Since  the  half-power  bandwidth  of  the  stop  band  is  1  GHz  at  0  = 
30°,  a  2-GHz  shift  in  the  resonant  frequency  over  the  range  of  incident  angles  on  the  array 
causes  serious  distortion  in  the  frequency  response.  Locally,  each  dipole  is  resonating  at  a 
different  frequency,  due  to  the  different  effective  incident  angles  at  each  point  on  the  array. 

As  a  remedy,  the  idea  is  to  design  each  cell  as  though  it  were  the  unit  cell  of  a  periodic 
surface  illuminated  by  a  plane  wave  at  the  effective  incident  angle  seen  by  that  particular 
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Figure  6.6.  Shift  in  resonant  frequency  for  the  element  of  Ott  et  al.  [37], 


cell.  In  the  15-element  example,  unit  cells  are  designed  for  15  periodic  surfaces  illuminated 
by  plane  waves  at  15  different  incident  angles,  such  that  they  all  resonate  at  the  same  center 
frequency.  The  15  unit  cells  obtained  in  this  manner  are  assembled  to  construct  the  15- 
element  TTay,  and  the  result  is  tested  in  a  finite  frequency  selective  surface.  In  the  design 
stage,  however,  only  the  horizontal  cell  sizes  and  element  lengths  are  permitted  to  vary, 
while  the  vertical  cell  size  Ty  and  element  width  are  held  constant. 

In  designing  the  15  cells  that  compose  the  finite  array,  two  methods  were  used.  A 
design  frequency  is  chosen,  at  which  all  of  the  elements  of  the  finite  array  were  to  become 
resonant.  In  the  first  method,  both  the  lengths  and  spacings  of  the  elements  were  scaled  by 
the  ratio  of  the  design  frequency  and  the  resonant  frequency  at  the  dipole's  effective 
illumination  angle.  In  the  second  method,  only  the  spacings  between  the  elements  were 
adjusted,  according  to  the  above  ratio,  leaving  the  lengths  fixed.  After  assembling  the  15 
unit  cells  obtained  in  this  manner  and  computing  their  frequency  response  in  the  finite 
array,  the  second  method  was  found  to  give  somewhat  better  performance  than  the  first 
method,  and  its  reflection  curve  is  shown  in  Fig.  6.7.  Hence,  as  a  result  of  tapering  the 
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surface  geometry,  all  of  the  power  is  reflected  at  resonance  as  before,  leaving  only  a  slight 
reduction  in  bandwidth,  which  may  be  recoverable  through  further  optimization. 


Nonuniform  Array 


Figure  6.7.  Power  reflection  coefficient  versus  frequency:  plane-wave 
scattering  from  a  periodic  surface,  and  scattering  from  a 
tapered  finite-width  FSS  of  15  elements  due  to  the 
electromagnetic-hom  incident  field. 
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CHAPTER  7 

CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  WORK 

The  analysis  of  finite  frequency  selective  surfaces  is  still  in  its  infancy:  this  thesis  is 
the  result  of  one  of  the  first  attempts  at  modeling  finite  arrays.  In  it,  the  formulation  of  the 
scattering  problem  in  the  spectral  domain  has  been  presented,  the  essential  considerations 
involved  in  implementing  the  spectral-Galerkin  equation  have  been  discussed,  and  a  few 
selected  geometries  have  been  analyzed  to  illustrate  the  conclusions.  Motivated  by  existing 
applications,  numerical  modeling  of  generalized  sources  and  nonuniform  frequency 
selective  surfaces  has  also  been  done. 

The  effects  of  finite  dimensions  on  the  performance  of  frequency  selective  surfaces 
have  been  classified  into  two  frequency  regimes.  At  frequencies  near  the  first  resonance  of 
the  element,  the  edge  effect  is  localized  and  extends  from  three  to  seven  cells  toward  the 
center  of  the  array.  Sample  ca-culations  show  that  the  component  of  the  current  flowing 
perpendicularly  to  the  polarization  of  the  incident  field  is  significantly  enhanced  near  the 
edges;  however,  for  more  than  1 5  elements,  the  periodic  approximation  gives  about  two 
per  cent  accuracy  in  the  scattered  fields.  At  lower  frequencies,  an  edge-to-edge  resonance 
can  be  excited  on  the  array,  even  when  there  is  no  direct  electrical  connection  between  the 
two  edges.  In  this  regime,  the  edge  effect  is  obviously  not  localized,  and  it  would  be 
impossible  to  accurately  predict  the  scattered  fields  from  such  a  surface  using  infinite 
models. 

In  Chapter  6,  it  was  shown  that  the  use  of  a  realistic  source  model,  such  as  the  near 
field  of  an  electromagnetic  horn,  caused  serious  degradation  in  the  frequency  response  of  a 
particular  finite  frequency  selective  surface.  Since  the  resonant  frequency  of  most 
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frequency  selective  surfaces  is  strongly  influenced  by  incident  angles,  it  was  postulated  that 
the  poor  performance  was  due  to  the  presence  of  a  wide  range  of  effective  incident  angles 
on  the  array.  A  method  was  therefore  suggested  to  continuously  vary  the  lattice  of  the 
array  in  order  to  force  all  of  the  elements  to  resonate  simultaneously  at  a  given  frequency. 
Based  on  this  theory,  a  nonuniform  frequency  selective  surface  was  designed  and  modeled 
in  the  presence  of  the  horn,  and  it  was  shown  that  continuously  varying  the  surface 
geometry  was  an  effective  means  of  alleviating  the  distortion. 

Looking  ahead  to  possible  areas  of  research  in  frequency  selective  surfaces,  it  is 
appropriate  to  mention  first  of  all  that  this  is  not  the  last  word  in  finite  frequency  selective 
surfaces.  Because  of  the  large  amounts  of  computer  time  required  for  some  geometries  and 
limitatiohs  in  computer  storage,  the  two-dimensional  examples  were  limited  to  certain  thin- 
patch  geometries  that  could  be  analyzed  with  a  small  number  of  unknowns.  The  next 
generation  of  computers  will  make  many  of  the  assumptions  used  here  unnecessary. 
Furthermore,  all  of  the  results  presented  in  the  thesis  were  generated  with  special-purpose 
routines.  A  general-purpose  program  capable  of  modeling  arbitrary  finite  frequency 
selective  surfaces  with  arbitrary  incident  fields  would  be  a  valuable  contribution. 

The  simplifications  and  restrictions  imposed  at  the  beginning  suggest  other  possible 
areas  of  research.  Many  applications  that  require  wider  frequency  bandwidths  will  need  to 
incorporate  two  or  more  printed  surfaces  to  meet  those  requirements.  While  this  problem 
has  been  solved  satisfactorily  for  periodic  surfaces,  numerical  modeling  of  multiple  layers 
of  finite  frequency  selective  surfaces  has  never  been  attempted.  Computation  time  and 
storage  requirements  increase  in  proportion  to  the  number  of  layers  squared,  if  the 
multilayered  problem  is  solved  at  once.  However,  a  procedure  analogous  to  the  scattering 
matrix  approach  for  periodic  surfaces  would  reduce  the  demands  on  computer  resources. 

A  second  area  of  research  is  suggested  by  the  fact  that  in  this  work,  only  the  printed 
surface  was  assumed  to  be  finite  in  some  sense,  while  the  surrounding  medium  was 
assumed  to  be  an  infinite  medium  of  uniform  homogeneous  dielectric  layers.  This  made  it 
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possible  to  derive  the  analytical  Green’s  functions.  If  the  guided  modes  in  the  dielectric 
wave  guide  carry  sufficient  energy,  panicularly  for  oblique  incidence,  then  truncating  the 
dielectric  may  have  a  significant  effect  on  the  scattered  fields.  How  much  effect,  or  under 
what  condiuons  this  effect  may  be  substantial,  is  not  known.  Incorporating  the  effects  of  a 
finite  dielectric  slab  will  involve  the  equivalent  of  discretizing  the  polarization  current  in  the 
dielectric  medium  and  solving  simultaneously  for  this  current  and  the  conduction  current  on 
the  printed  surface.  Some  of  this  work  has  already  begun  at  the  University  of  Illinois  by 
considenng  a  finite  number  of  infinite  strips  embedded  in  an  infinite  dielectric  cylinder  of 
rectangular  cross  section. 

A  third  area  of  research  is  that  of  curved  surfaces.  A  natural  starting  place  would  be 
the  development  of  a  formal  notation  capable  of  describing  the  shape,  position,  orientation, 
and  curvature  of  an  element  in  three-dimensional  space,  much  like  the  notation  used  in  the 
geometrical  theory  of  diffraction.  Such  a  formalism  might  also  be  applied  to  other 
deterministic  scattering  problems  as  well. 

Finally,  a  fourth  area  of  research  that  could  be  carried  out  in  the  context  of  periodic 
surfaces  is  the  development  of  a  systematic  design  procedure.  Most,  if  not  all,  of  the 
programs  used  in  practice  are  analysis  tools.  The  designer  is  given  a  set  of  specifications, 
and  then,  based  on  his  experience,  he  analyzes  tens  or  hundreds  of  cases  to  find  one  that 
comes  close  to  meeting  the  specifications.  Depending  on  his  level  of  expertise,  this  process 
could  be  very  time-consuming  -  and  the  intuition  he  gains  is  not  easily  passed  on  to  his 
colleagues. 

The  alternative  would  be  to  automate  this  procedure  in  an  optimization  program.  TTie 
inputs  to  the  program  would  be  the  design  specifications;  resonant  frequency,  bandwidth, 
polarization  of  the  field  (linear  or  circular),  etc.  Beginning  at  an  initial  design,  the  program 
would  be  designed  to  incrementally  vary  the  adjustable  parameters  of  the  surface  until  the 
design  goals  were  achieved.  The  output  of  the  program  would  be  the  number  of 
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conducting  layers,  the  lattice,  the  shapes  and  sizes  of  the  elements,  and  the  thicknesses  and 
permittivities  of  the  dielectric  layers. 

Predicting  the  future  is  a  tricky  business,  but  the  problems  discussed  here  are 
potentially  interesting  and  are  motivated  by  real  applications  in  existing  or  proposed 
systems. 
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